//  only for private communication
//
//  Xianguo Lu 
//  lu@physi.uni-heidelberg.de
//  Xianguo.Lu@cern.ch
//

#include <fstream>

#include "Math/Functor.h"
#include "Math/Factory.h"
#include "Math/Minimizer.h"

#include "TASImage.h"
#include "TCanvas.h"
#include "TColor.h"
#include "TChain.h"
#include "TF1.h"
#include "TFile.h"
#include "TGaxis.h"
#include "TGraph.h"
#include "TGraphAsymmErrors.h"
#include "TGraphErrors.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TLine.h"
#include "TMath.h"
#include "TMatrixD.h"
#include "TMinuit.h"
#include "TLatex.h"
#include "TLegend.h"
#include "TPaletteAxis.h"
#include "TRandom3.h"
#include "TROOT.h"
#include "TStopwatch.h"
#include "TStyle.h"
#include "TTree.h"
#include "TVector3.h"

#include "TTreeStream.h"

#include "XGLUtils.h"

#define EPSILON 1e-12

Int_t XGLUtils::fgFitPrintLevel = 1;
Int_t XGLUtils::fgFitNStep=0;
Int_t XGLUtils::fgFitNData=0;
Double_t * XGLUtils::fgFitX = 0x0;
Double_t * XGLUtils::fgFitY = 0x0;
Double_t * XGLUtils::fgFitEx = 0x0;
Double_t * XGLUtils::fgFitEy = 0x0;
XGLUtils::FFunc XGLUtils::fgFitFunc = 0x0;
Double_t XGLUtils::fTol = 1;
Int_t XGLUtils::fMaxIter = (Int_t)1e8;
Int_t XGLUtils::fMinOpt = XGLUtils::kTMINUIT;

//===================================================================  
//                               Algorithm 
//===================================================================
Double_t XGLUtils::GetCorrelatedError(const Double_t v0, const Double_t v1, const Double_t a0, const Double_t a1, const Double_t a2)
{
  //                   ( a0   a1 ) (v0)
  //df^2 = (v0,v1) *   ( a1   a2 ) (v1) = v0*v0*a0 + v1*v1*a2 + 2*v0*v1*a1
  //f=f(a,b), v0=df/da, v1=df/db
  //f=a/b, v0=f/a, v1=-f/b

  const Double_t df2 = v0*v0*a0 + v1*v1*a2 + 2*v0*v1*a1;
  /*
  //XGLUtils::GetCorrelatedError df2<0 -2.980232e-08 1.000000e-10 -- 1.640908e+07 -1.640908e+07 -- 3.221823e-07 3.221823e-07 3.221823e-07
  const Double_t eps = 1e-2;
  if(df2<eps){
    if(TMath::Abs(df2)<eps){
      return 0;
    }
    else{
      printf("XGLUtils::GetCorrelatedError df2<0 %e %e -- %e %e -- %e %e %e\n", df2, eps, v0, v1, a0, a1, a2);exit(1);
    }
  }
  */
  if(df2<EPSILON){
    return 0;
  }
  return TMath::Sqrt( df2 );
}

Double_t XGLUtils::GetMeanRMS(const Double_t nn, const Double_t sum, const Double_t w2s, Double_t * grms, Double_t * gerr)
{
  //
  //calculate mean (with error) and rms from sum, w2s, nn
  //if nn=0, mean, error, and rms all = 0
  //

  Double_t tmean = 0, trms = 0, terr = 0;

  const Double_t eps = 1e-10;
  if(nn>eps){
    tmean = sum/nn;

    const Double_t arg = w2s/nn-tmean*tmean;
    if(TMath::Abs(arg)<eps){
      trms = 0;
    }
    else{
      if( arg <0 ){
        printf("XGLUtils::GetMeanRMS error negative sqrt argument!! %e -- %e %e %f\n", arg, w2s, sum, nn); exit(1);
      }
      
      trms = TMath::Sqrt(arg);
    }

    terr = trms/TMath::Sqrt(nn);
  }

  if(grms){
    (*grms) = trms;
  }

  if(gerr){
    (*gerr) = terr;
  }

  return tmean;

}

Double_t XGLUtils::TruncatedMean(TTree *tt, const TString var, const Double_t lowfrac, const Double_t highfrac, Double_t * grms, Double_t * gerr)
{
  const Int_t nx = tt->GetEntries();
  Double_t xdata[nx];

  Double_t xx = -999;
  tt->SetBranchAddress(var, &xx);

  for(Int_t ii=0; ii<nx; ii++){
    tt->GetEntry(ii);
    xdata[ii]=xx;
  }

  return TruncatedMean(nx, xdata, lowfrac, highfrac, grms, gerr);
}

Double_t XGLUtils::TruncatedMean(const Int_t nx, const Double_t xdata[], const Double_t lowfrac, const Double_t highfrac, Double_t * grms, Double_t * gerr, Double_t *wws)
{
  //
  //calculate truncated mean
  //return <x*w>_{low-high according to x}
  //

  /*
  //test->
  for(Int_t ii=0; ii<nx; ii++){
    printf("test %d/%d %f\n", ii, nx, xdata[ii]);
  }
  //<--test
  */

  Int_t index[nx];
  TMath::Sort(nx, xdata, index, kFALSE);

  Int_t nused = 0;
  Double_t sum = 0;
  Double_t w2s = 0;
  const Int_t istart = Int_t (nx*lowfrac);
  const Int_t istop  = Int_t (nx*highfrac);

  //=,< correct, because when low=0, high=1 it is correct
  for(Int_t ii=istart; ii<istop; ii++){
    Double_t weight = 1;
    if(wws){
      weight = wws[index[ii]];
    }
    const Double_t sx = xdata[index[ii]]*weight;

    sum += sx;
    w2s += sx*sx;

    nused++;
    //printf("test in loop %d/%d %f %f %f\n", ii, nused, sx, sum, w2s);
    
  }

  return GetMeanRMS(nused, sum, w2s, grms, gerr);

}

//Numerical Recipes in C, 1997, 3.1 Polynomial Interpolation and Extrapolation, p109
//Neville's algorithm
//error defined here depends on the choice of the route from ya[ns] to y_final
void XGLUtils::Polint(const Int_t n, const Double_t *xa, const Double_t *ya, const Double_t x, Double_t &y, Double_t &dy, const Int_t kprint)
{
  /*
  for(Int_t ii=0; ii<n; ii++){
    printf("XGLUtils::Polint %10.10e %03d/%03d: %10.10e %10.10e\n", x, ii, n, xa[ii], ya[ii]);
  }
  */

  Double_t c[n];
  Double_t d[n];
  Double_t dift=-999;

  Int_t ns=0;
  Double_t dif = fabs(x-xa[ns]); 
  for (Int_t i=0; i<n;i++) {
    if ( (dift=fabs(x-xa[i])) < dif) { 
      ns=i; 
      dif=dift; 
    }
    c[i]=ya[i]; 
    d[i]=ya[i];

    if(kprint) printf("polint =====> x[%d]: %3.10e\ty[%d]: %3.10e\n", i, xa[i], i, ya[i]);
  }

  y = ya[ns];
  if(kprint) printf("======> y0: %3.10e at ns=%d\n", y, ns);

  Int_t did=ns-1;
  for (Int_t m=1;m<=n-1;m++){
    for (Int_t i=0; i<n-m; i++){
      Double_t ho=xa[i]-x; 
      Double_t hp=xa[i+m] -x;
      Double_t w = c[i+1]-d[i];
      Double_t den=ho-hp;
      if ( den == 0.0){
        printf("XGLUtils::Polint error (den=ho-hp) == 0.0!! %e %e %e\n", x, xa[i], xa[i+m]); exit(1);
      }

      c[i]=ho*w/den;
      d[i]=hp*w/den; 
      //if(kprint) printf("======> m %d, c[%d]: %5.6e\t d[%d]: %5.6e\n", m,i,c[i],i, d[i]);
    }
    if( m <= n -1-ns ){
      dy = c[ns];
      if(kprint) printf("------------> m %d, c[%d]: %5.6e\n", m, ns, dy);
    }
    else {
      dy = d[did];
      if(kprint) printf("------------> m %d, d[%d]: %5.6e\n", m, did, dy);
      did--;
    }
    y += dy;
    if(kprint) printf("======> y%d: %3.10e\n",m, y);
  }
}

Bool_t XGLUtils::FindY(FFunc func, Double_t *par, const Double_t xmin, const Double_t xmax, const Double_t ytarget, Double_t &pos, Int_t &ndepth, const Double_t eps, const Double_t resolution)
{
  ndepth++;

  Bool_t kFindMax = kFALSE;
  Bool_t kFindMin = kFALSE;
  Bool_t kFindY   = kFALSE;
  if(ytarget>1e9)
    kFindMax = kTRUE;
  else if(ytarget<-1e9)
    kFindMin = kTRUE;
  else
    kFindY = kTRUE;

  const Double_t old = pos;

  Double_t x0 = xmin;
  Double_t dist =  1e10;
  Double_t ymax = -1e10;
  Double_t ymin =  1e10;
  Double_t b0 = -999, b1 = -999;
  while(x0<xmax){
    Double_t x1 = x0 + resolution;
    if(x1>xmax)
      x1 = xmax;

    const Double_t xevl = x0 + 0.5*(x1-x0);
    const Double_t yevl = func(&xevl, par);
    //printf("jjjjj %.10e %.10e %.10e %.20e\n", x0, x1, xevl, yevl);
    
    Bool_t kUpdateMax = kFALSE;
    Bool_t kUpdateMin = kFALSE;
    Bool_t kUpdateDist = kFALSE;
    if(yevl>ymax){
      kUpdateMax = kTRUE;
      ymax = yevl;
    }
    if(yevl<ymin){
      kUpdateMin = kTRUE;
      ymin = yevl;
    }
    if( TMath::Abs(yevl-ytarget)<dist ){
      kUpdateDist = kTRUE;
      dist = TMath::Abs(yevl-ytarget);
    }

    if( 
       (kFindMax && kUpdateMax)
       || (kFindMin && kUpdateMin)
       || (kFindY && kUpdateDist)
        ){
      pos = xevl;
      b0 = x0;
      b1 = x1;
    }
    x0 += resolution;
  }

  if(TMath::Abs(old-pos)>eps){
    if(ndepth>=50){
      printf("XGLUtils::FindY ndepth >= 50!! %d: %.10e %.10e %.10e, %.10e, %e %.10e -- %e %e %e\n", ndepth, ytarget, xmin, xmax, pos, eps, resolution, old, pos, TMath::Abs(old-pos));
      return kFALSE;
    }
    else
      return FindY(func, par, b0, b1, ytarget, pos, ndepth, eps, pow(resolution,2)/(xmax-xmin));
  }
  return kTRUE;
}

Bool_t XGLUtils::GetMPVFWHM(FFunc func, Double_t *par, const Double_t xmin, const Double_t xmax, Double_t *vout, const Double_t eps, Double_t resolution)
{
  if(resolution<0)
    resolution = (xmax-xmin)/100;

  Double_t ytarget = 1e10; //>1e9 for maximum
  Double_t pos = -999;
  Int_t ndepth = 0;
  if(!FindY(func, par, xmin, xmax, ytarget, pos, ndepth, eps, resolution)){
    printf("XGLUtils::GetMPVFWHM FindY(1) fail! %f %f %e %e %d %e %e\n", xmin, xmax, ytarget, pos, ndepth, eps, resolution);
    return kFALSE;
  }
  const Double_t mpv = pos;

  //----------------- 

  const Double_t ymax = func(&mpv, par);
  ytarget = ymax/2;

  pos = -999;
  ndepth = 0;
  if(!FindY(func, par, xmin, mpv, ytarget, pos, ndepth, eps, resolution)){
    printf("XGLUtils::GetMPVFWHM FindY(2) fail! %f %f %e %e %d %e %e\n", xmin, mpv, ytarget, pos, ndepth, eps, resolution);
    return kFALSE;
  }

  const Double_t leg0 = pos;

  pos = -999;
  ndepth = 0;
  if(!FindY(func, par, mpv, xmax, ytarget, pos, ndepth, eps, resolution)){
    printf("XGLUtils::GetMPVFWHM FindY(3) fail! %f %f %e %e %d %e %e\n", mpv, xmax, ytarget, pos, ndepth, eps, resolution);
    return kFALSE;
  }

  const Double_t leg1 = pos;
  const Double_t fwhm = leg1-leg0;
  vout[0]=mpv;
  vout[1]=fwhm;
  vout[2]=leg0;
  vout[3]=leg1;

  return kTRUE;
}

//===================================================================
//                                Functions
//===================================================================
Double_t XGLUtils::ProtectedExp(const Double_t xx)
{
  /*
    root [14] exp(709)
    (const double)8.21840746155497238e+307
    root [15] exp(710)
    Error: exp param[0]=710 up:709 low:-inf out of range (tmpfile):1:
  */
  const Double_t maxexparg = 700;
  
  return TMath::Exp(TMath::Min(maxexparg, xx));
}

Double_t XGLUtils::ProtectedPow(const Double_t xx, const Double_t yy)
{
  if(xx<0){
    printf("XGLUtils::ProtectedPow xx<0 %e\n", xx); exit(1);
  }
  //xx^yy = exp(log(xx^yy))=exp(yy log(xx))
  return ProtectedExp(yy*TMath::Log(xx+EPSILON));
}

Double_t XGLUtils::QuadSum(const Double_t x1, const Double_t x2)
{
  return TMath::Sqrt( TMath::Power(x1,2) + TMath::Power(x2,2) );
}

Double_t XGLUtils::QuadSum(const Double_t x1, const Double_t x2, const Double_t x3)
{
  return TMath::Sqrt( TMath::Power(x1,2) + TMath::Power(x2,2) + TMath::Power(x3,2) );
}

Double_t XGLUtils::ALG(const Double_t * xx, const Double_t * par)
{
  //
  //npar = 4
  //
  const Double_t eps = 1e-6;

  const Int_t ng = 5;
  const Double_t xmin = TMath::Max(0.0, xx[0]-ng*par[3]);
  const Double_t xmax = xx[0]+ng*par[3];

  TF1 ftmp("xgltuilsalg",Form("[0]*TMath::Landau(x, [1], [2], 1)*TMath::Gaus(x, %.10f, [3], 1)", xx[0]), xmin, xmax);

  const Double_t sum = ftmp.Integral(xmin, xmax, par, eps);
  
  return sum;
}

Double_t XGLUtils::NorDimlessLGE(const Double_t * par)
{
  const Double_t xmin = 0;
  const Double_t xmax = 100;

  const Double_t eps = 1e-6;

  TF1 ftmp1("xglutilsnordimndl1","TMath::Landau(  x, [0], [1]*[0], 1)*TMath::Exp(-[3]*    x/[0])",     xmin,   xmax);
  TF1 ftmp2("xglutilsnordimndl2","TMath::Landau(1/x, [0], [1]*[0], 1)*TMath::Exp(-[3]*(1/x)/[0])/x/x",    0, 1/xmax);

  const Double_t nor1 = ftmp1.Integral(xmin, xmax, par, eps);
  const Double_t nor2 = ftmp2.Integral(0, 1/xmax, par, eps);
  return nor1+nor2;
}

Double_t XGLUtils::DimlessALGE(const Double_t * xx, Double_t * par)
{
  const Int_t npar = 5;

  const Double_t eps = 1e-6;
  Double_t absp[npar];
  for(Int_t ii=0; ii<npar; ii++)
    absp[ii] = TMath::Abs(par[ii]);

  const Double_t gs = absp[3]*absp[1];
  const Int_t ng = 5;
  const Double_t xmin = TMath::Max(0.0, xx[0]-ng*gs);
  const Double_t xmax = xx[0]+ng*gs;

  TF1 ftmp("xglutilsdimlessalgefunc",Form("[0]*TMath::Landau(x, [1], [2]*[1], 1)*TMath::Gaus(x, %.10f, [3]*[1], 1)*TMath::Exp(-[4]*x/[1])", xx[0]), xmin, xmax);

  const Double_t sum = ftmp.Integral(xmin, xmax, absp, eps);
  const Double_t nor = NorDimlessLGE(&(absp[1]));

  return sum/(nor+EPSILON);
}

Double_t XGLUtils::InvSqrt(const Double_t *xx, const Double_t *par)
{
  if(xx[0]<EPSILON){
    printf("XGLUtils::InvSqrt x<epsilon!!\n"); exit(1);
  }
  return par[0]/TMath::Sqrt(xx[0]);
}

Double_t XGLUtils::PolN(const Double_t *xx, const Double_t *par, const Int_t nn)
{
  //
  //
  //
  Double_t sum = 0;
  for(Int_t ii=0; ii<=nn; ii++){
    sum += par[ii]*TMath::Power(xx[0],ii);
  }
  return sum;
}
Double_t XGLUtils::Pol1(const Double_t *xx, const Double_t *par)
{
  return PolN(xx, par, 1); 
}

Double_t XGLUtils::Pol2(const Double_t *xx, const Double_t *par)
{
  return PolN(xx, par, 2);
}

Double_t XGLUtils::Pol3(const Double_t *xx, const Double_t *par)
{
  return PolN(xx, par, 3);
}

Double_t XGLUtils::ALandau(const Double_t *xx, const Double_t *par)
{
  //
  //p0*landau, npar = 3
  //
  return par[0]*Landau(xx, &(par[1]));
}

/*
Double_t XGLUtils::NorExpLandau(Double_t *par)
{
  //
  //normalization of explandau
  //

  const Double_t xmin = 0;
  const Double_t xmax = 100;

  TF1 ftmp1("tmpexplandaundl1","TMath::Landau(  x, [0], TMath::Abs([1]), 1)*TMath::Exp(-TMath::Abs([2])*    x)",     xmin,   xmax);
  TF1 ftmp2("tmpexplandaundl2","TMath::Landau(1/x, [0], TMath::Abs([1]), 1)*TMath::Exp(-TMath::Abs([2])*(1/x))/x/x",    0, 1/xmax);

  //1e-6 best 
  const Double_t eps = 1e-3;
  const Double_t nor1 = ftmp1.Integral(xmin, xmax, par, eps);
  const Double_t nor2 = ftmp2.Integral(0, 1/xmax, par, eps);

  const Double_t totn = nor1+nor2;

  return totn;
}

Double_t XGLUtils::ExpLandau(const Double_t *xx, Double_t *par)
{
  //
  //exp*landau
  //
  const Double_t totn = NorExpLandau(par);
  return TMath::Exp(-TMath::Abs(par[2])*xx[0])*Landau(xx, par)/(totn+EPSILON);
}

Double_t XGLUtils::DimLessExpLandau(const Double_t *xx, const Double_t *pp)
{
  //
  //exp*landau, mpv = p0, wid = p1*p0, kappa = p2/p0
  //
  const Double_t par[]={pp[0], pp[1]*pp[0], pp[2]/(pp[0]+EPSILON)};

  return ExpLandau(xx, par);
}
*/

Double_t XGLUtils::InvLandau(const Double_t *yy, const Double_t *par)
{
  //
  //Landau(x)dx -> f(y) dy, y=1/x, f(y(x))=Landau(x)*x*x
  //f(y)= Landau(1/y)/y/y;
  //

  const Double_t xx[1]={1./(yy[0]+EPSILON)};
  return Landau(xx, par)/yy[0]/yy[0];
}

Double_t XGLUtils::Landau(const Double_t *xx, const Double_t *par)
{
  //
  //landau, npar = 2
  //
  return TMath::Landau(xx[0], par[0], TMath::Abs(par[1]), kTRUE);
}

Double_t XGLUtils::ResolutionGaus(const Double_t *xx, const Double_t *par)
{
  //npar = 3

  //
  //root [0] TMath::Gaus(3, 0, -9,1)
  //(Double_t)(-4.19314697436659056e-02)
  //root [1] TMath::Gaus(3, 0, 9,1)
  //(Double_t)4.19314697436659056e-02
  //
  //!! width has to be abs-ed!!
  return par[0]*TMath::Gaus(xx[0], par[1], par[1]*TMath::Abs(par[2]), kTRUE);
}

Double_t XGLUtils::RawGaus(const Double_t *xx, const Double_t *par)
{
  //npar = 3

  return par[0]*TMath::Gaus(xx[0], par[1], TMath::Abs(par[2]), kTRUE);
}

/*
Double_t XGLUtils::ExpGaus(const Double_t xx[], const Double_t par[])
{
  const Double_t amp    = TMath::Abs(par[0]);
  const Double_t mean   = TMath::Abs(par[1]);
  const Double_t sigma  = TMath::Abs(par[2]);
  const Double_t lambda = TMath::Abs(par[3]);

  
  // //http://www.wolframalpha.com/input/?_=1348047374797&fp=1&i=integrate+exp%28-%28x-m%29^2%2fs^2%2f2+-+x+L%29%2c+x+from+-inf+to+inf&incTime=true
  // const Double_t normalization = TMath::Sqrt(2 * TMath::Pi()) * sigma * TMath::Exp(0.5* lambda *( lambda * sigma*sigma -  2 * mean));

  // if(normalization>EPSILON)
  //   return amp*TMath::Exp(-TMath::Power((xx[0]-mean)/sigma,2)/2. - xx[0]*lambda)/normalization;
  

  //http://www.wolframalpha.com/input/?_=1348048413665&i=integrate+exp%28-%28x-m%29^2%2fs^2%2f2+-+%28x-m%29+L%29%2c+x+from+-inf+to+inf&fp=1&incTime=true  
  return amp*TMath::Exp(-TMath::Power((xx[0]-mean)/sigma,2)/2. - (xx[0]-mean)*lambda - 0.5* lambda * lambda * sigma*sigma)/ TMath::Sqrt(2 * TMath::Pi()) * sigma;
}
*/
Double_t XGLUtils::TailedGaus0(const Double_t xx, const Double_t amp, const Double_t mean, const Double_t sigma, const Double_t lambda)
{
  //
  //exp(-x/lambda)
  //
  
  Double_t sum=0, sW=0;

  //lambda=0, exp(-x/lambda)/lambda -> delta function
  if(lambda<EPSILON){
    return amp*TMath::Gaus( xx, mean, sigma, kTRUE);
  }

  //nlam < 5 step not enough to cover to tail to 1e-4
  //nlam < 10 not enough to 1e-5
  for(Double_t nlam=0; nlam<15; nlam+=0.5){
    const Double_t weight = TMath::Exp(-nlam);
    sum += TMath::Gaus( xx - nlam*lambda, mean, sigma, kTRUE)*weight;
    sW+=weight;
  }

  /*
  //exp(-x/lambda)
  const Double_t step = 0.1;
  Double_t weight = 1;
  Double_t nlam = 0;
  const Double_t eps = 1e-9;
  while(weight>eps){
    weight = TMath::Exp(-nlam);
    const Double_t delta = TMath::Gaus( xx - nlam*lambda, mean, sigma, kTRUE)*weight;
    sum += delta;
    sW += weight;
    nlam+=step;
  }
  */

  /*
  const Double_t kappa = 1/lambda;
  const Double_t eps = 1e-3;
  const Double_t nstep = 20;
  Double_t tau = 0;
  for(Int_t ii=0; ii<nstep; ii++){
    sum+= TMath::Gaus(xx-tau, mean, sigma,kTRUE);
    
    //const Double_t oldx = tau;
    tau = -TMath::Log(TMath::Exp(-kappa*tau)-(1-eps/kappa)/nstep)/kappa;
    //printf("%d %e %e -- %e -- %e\n", ii,oldx, tau, exp(-kappa*tau)*kappa, exp(-kappa*oldx)*kappa-exp(-kappa*tau)*kappa);
    sW += 1;
  }
  */
  /*
  const Double_t eps = 1e-3;
  const Double_t nstep = 20;
  //exp(-id*sigma/nstep/lambda)<eps
  const Double_t NN = nstep*sigma/lambda/TMath::Log(1/eps);

  Double_t weight = 1;

  const Double_t dx = sigma/100;
  Double_t id = 0;
  while(weight>eps){
    const Double_t tau = id*dx;
    weight = TMath::Exp(-tau/lambda);
    sum += TMath::Gaus( xx - tau, mean, sigma, kTRUE)*weight;
    sW+=weight;
    id++;
  }
  */

  /*
  sW = 1;
  const Double_t kappa = 1/lambda;
  const Double_t eps = 1e-6;
  const Double_t nstep = 20;
  Double_t tau = 0;
  for(Int_t ii=1; ii<=nstep; ii++){
    //const Double_t oldx = tau;
    tau = -TMath::Log(TMath::Exp(-kappa*tau)-(1-eps/kappa)/nstep)/kappa;
    sum+= TMath::Gaus(xx-tau, mean, sigma,kTRUE);
    //printf("%d %e %e -- %e -- %e\n", ii,oldx, tau, exp(-kappa*tau)*kappa, exp(-kappa*oldx)*kappa-exp(-kappa*tau)*kappa);
  }
  */

  if (sW>EPSILON) return amp*sum/sW;
  return 0;
}

Double_t XGLUtils::TailedGaus(const Double_t xx[], const Double_t par[])
{
  const Double_t amp    = TMath::Abs(par[0]);
  const Double_t mean   = TMath::Abs(par[1]);
  const Double_t resolution  = TMath::Abs(par[2]);
  const Double_t lambda = TMath::Abs(par[3]);

  return TailedGaus0(xx[0], amp, mean, resolution*mean, lambda);
}

Double_t XGLUtils::GausN(const Double_t *xx, const Double_t *par, const Int_t ntype, FFunc basefunc, const Int_t npartype)
{
  //npar = ntype*npartype

  Double_t sum = 0;

  for(Int_t itype=0; itype<ntype; itype++){
    sum += basefunc(xx, &(par[npartype*itype]));
  }

  return sum;
}

Double_t XGLUtils::HadronBetheBloch(const Double_t *xx, const Double_t *par)
{
  //RawALEPH has numerical difficulty in ln(p2+1/pow(bg, p4))
  //use flat dedx for electron. deviation from this at high bg is not inportant because anyway there the statistics is very limited
  //npar = 5;

  //hadronBB npar = 4

  const Double_t bg = xx[0];
  const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);

  const Double_t p0 = par[0];//overall prefactor
  const Double_t p1 = par[1];//beta exponent ~ 2
  const Double_t p2 = par[2];
  const Double_t p3 = par[3];

  const Double_t bethebloch = p0 *( 1 + p2*TMath::Log(bg) + p3*TMath::Log(beta) )/TMath::Power(beta, p1);

  return bethebloch;
}

Double_t XGLUtils::FullBetheBloch(const Double_t *xx, const Double_t *par)
{
  //RawALEPH has numerical difficulty in ln(p2+1/pow(bg, p4))
  //use flat dedx for electron. deviation from this at high bg is not inportant because anyway there the statistics is very limited
  //npar = 5;

  //fullBB npar = 5

  const Double_t bg = xx[0];
  const Double_t p4 = par[4];//plateau ~ 80

  const Double_t bethebloch = HadronBetheBloch(xx, par);

  const Double_t tp = 10;
  if(bg<tp){
    return bethebloch;
  }
  else{
    return TMath::Min(p4, bethebloch);
  }
}

Double_t XGLUtils::RawAleph(const Double_t *xx, const Double_t *par)
{
  const Double_t bg = xx[0];
  const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);

  const Double_t p0 = TMath::Abs(par[0]);
  const Double_t p1 = TMath::Abs(par[1]);
  const Double_t p2 = TMath::Abs(par[2]);
  const Double_t p3 = TMath::Abs(par[3]);
  const Double_t p4 = TMath::Abs(par[4]);

  //return p0*((p1-TMath::Log(p2+1/TMath::Power(bg,p4)))/TMath::Power(beta,p3) - 1 );
  return p0*( (p1- TMath::Log(p2)- TMath::Log(1+1/(TMath::Power(bg,p4)*p2))  )/TMath::Power(beta,p3) - 1 );
}

Double_t XGLUtils::Lund(const Double_t * xx,  const Double_t * rawpar)
{
  Double_t par[100];
  memcpy(&(par[1]), rawpar, sizeof(Double_t)*5);

  const Double_t bg = xx[0];

  // bg is beta*gamma
  const Double_t beta2 = bg*bg / (1.0+bg*bg);
  const Double_t a = par[1];
  const Double_t b = par[2]*par[1];//in this way value = rawpar[0]*()
  const Double_t c = par[3];
  const Double_t e = par[4];
  const Double_t f = par[5];
  const Double_t d = TMath::Exp(c*(a-f)/b);
  
  const Double_t powbg = TMath::Power(1.0+bg, c);
  
  const Double_t value = a/TMath::Power(beta2,e) +
    b/c*TMath::Log(powbg/(1.0 + d*powbg));
  return value;
}

//===================================================================
//                                Fitting 
//===================================================================
TMatrixDSym * XGLUtils::GetCov(TMinuit * tmn, const Bool_t ifixed[], const Double_t errs[])
{
  // case of fixed params need to take care
  //http://root.cern.ch/root/html/src/TMinuitMinimizer.cxx.html#Ilx.Z

  const Int_t nfree = tmn->GetNumFreePars();
  const Int_t npar = tmn->GetNumPars();

  TMatrixDSym matfree(nfree); 
  tmn->mnemat(matfree.GetMatrixArray(), nfree);

  TMatrixDSym *matall = new TMatrixDSym(npar);

  //r: row, c: column
  Int_t rfree = 0;
  for(Int_t rr=0; rr<npar; rr++){
    if(ifixed[rr]){
      continue;
    }
    Int_t cfree = 0;
    for(Int_t cc=0; cc<npar; cc++){
      if(ifixed[cc]){
        continue;
      }
      (*matall)[rr][cc] = matfree[rfree][cfree];
      cfree++;
    }
    rfree++;
  }

  /*
  //test->
    for(int ii=0; ii<NparBin(); ii++){
    const Int_t id = NparBuffer()+ii;
    printf("ii %d: err %e %e -- %e\n", ii, fLowErrs[id], sqrt((*matall)[id][id]), fLowErrs[id]-sqrt((*matall)[id][id]));
    for(int jj=0; jj<NparBin(); jj++){
    const Int_t i2 = NparBuffer()+jj;
    printf("%e ", (*matall)[id][i2]);
    }
    printf("\n");
    }
    //<-test
  */

  //check consistency!
  const Double_t eps = 1e-1;
  for(Int_t ipar=0; ipar<npar; ipar++){
    const Double_t cov = (*matall)[ipar][ipar];
    const Double_t sigma = TMath::Sqrt(cov);
    if(TMath::Abs(errs[ipar]-sigma)/errs[ipar]>eps){
      printf("XGLUtils::GetCov warning!! --> err-sigma > eps %d: %e %e %e -- %e\n", ipar, errs[ipar], cov, sigma, eps);
    }
  }

  return matall;
}

Int_t XGLUtils::FitKernel(const Int_t kml,const Int_t npar, Double_t pars[], Double_t errs[], Double_t chi[], const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs)
{
  //
  //kml: 0: chi2, 1: maximum likelihood, 2: binned likelihood
  //
  static Int_t counter = 0;
  if(counter==0){
    printf("\n\n**************** XGLUtils::FitKernel fMinOpt %d **************** \n", fMinOpt);
  }

  counter++;

  if(fMinOpt==kTMINUIT){
    return TMinuitFitKernel(kml, npar, pars, errs, chi, pfix, cov, lows, highs);
  }
  else{
    printf("XGLUtils::FitKernel MinimizerKernel no long term support! not turned on!\n"); exit(1);

    std::string minName;
    if(fMinOpt==kMINUIT){
      minName="Minuit";//this is exactly the same as TMinuitFitKernel
    }
    else if(fMinOpt==kMINUIT2){
      minName="Minuit2";
    }
    else{
      printf("XGLUtils::FitKernel fMinOpt error!! %d\n", fMinOpt); exit(1);
    }
    return MinimizerFitKernel(minName, kml, npar, pars, errs, chi, pfix, cov, lows, highs);
  }
}

Int_t XGLUtils::MinimizerFitKernel(const std::string minName, const Int_t kml,const Int_t npar, Double_t pars[], Double_t errs[], Double_t chi[], const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs)
{
  TStopwatch stwa;
  if(fgFitPrintLevel>=0){
    stwa.Start(kTRUE);
  }

  std::string algoName;
  algoName="Migrad";

  ROOT::Math::Minimizer *minimizer = ROOT::Math::Factory::CreateMinimizer(minName, algoName); 
  
  Double_t (*fcnptr)(const Double_t*) = 0x0;
  switch(kml){
  case 0:
    fcnptr = &GetFCNChisquare;
    break;
  case 1:
    fcnptr = &GetFCNMaximumLikelihood;
    break;
  case 2:
    fcnptr = &GetFCNBinnedLikelihood;
    break;
  default:
    printf("XGLUtils::MinimizerFitKernel bad kml %d\n", kml);exit(1);
  }

  ROOT::Math::Functor func(fcnptr, npar);
  minimizer->SetFunction(func); 
  
  minimizer->SetPrintLevel(fgFitPrintLevel); 
  
  const Double_t errdef = kml? 0.5 : 1;
  minimizer->SetErrorDef(errdef);

  minimizer->SetStrategy(1);

  minimizer->SetMaxFunctionCalls(fMaxIter);
  minimizer->SetMaxIterations(fMaxIter);

  const Double_t tol = (minName=="Minuit2"?0.5:1)*fTol;//0.5 due to internal definition of edm in Minuit2
  minimizer->SetTolerance(tol);
  
  Int_t nfree=0;
  for(Int_t ipar=0;ipar<npar;ipar++){
    const Double_t start=pars[ipar];
    const Double_t ss=TMath::Abs(pars[ipar])>EPSILON ? TMath::Abs(pars[ipar])*1e-1 : 1e-3;
    const Double_t low= (lows && lows[ipar]>-1e9)?lows[ipar]:0;
    const Double_t hi= (highs && highs[ipar]>-1e9)?highs[ipar]:0;
 
    if(pfix && pfix[ipar]){
      minimizer->SetFixedVariable(ipar, Form("p%d",ipar), start);
    }
    else{
      if(low<hi){
        minimizer->SetLimitedVariable(ipar, Form("p%d",ipar), start, ss, low, hi);
      }
      else{
        minimizer->SetVariable(ipar, Form("p%d",ipar), start, ss);
      }
      nfree++;
    }
  }
  
  //Minuit2-Migrad
  //http://root.cern.ch/root/html/ROOT__Minuit2__Minuit2Minimizer.html
  //status = 1    : Covariance was made pos defined
  //status = 2    : Hesse is invalid
  //status = 3    : Edm is above max
  //status = 4    : Reached call limit
  //status = 5    : Any other failure

  minimizer->Minimize(); 
  Int_t migradflag = minimizer->Status();
  
  //get par right after MIGRAD before HESSE modifies them
  Double_t migp[10];
  memcpy(migp, minimizer->X(), sizeof(Double_t)*npar);

  minimizer->Hesse();
  const Int_t hesseflag = minimizer->Status()-migradflag;

  //only correct migradflag after hesseflag set
  if(minName=="Minuit2" && migradflag==1){
    //pos def not required here
    migradflag=0;
  }

  const Int_t kfail = migradflag+hesseflag;
  if(kfail){
    printf("XGLUtils::FitKernel fitting error !! %d %d keep raw par and set err to -999!!!\n", migradflag, hesseflag);
    if(errs){
      for(Int_t ipar=0;ipar<npar;ipar++){
        errs[ipar]=-999;
      }
    }
    if(chi){
      chi[0]=chi[1]=-999;
    }
  }
  else{
    if(errs){
      memcpy(errs, minimizer->Errors(), sizeof(Double_t)*npar);
    }
    memcpy(pars, migp, sizeof(Double_t)*npar);
    if(chi){
      chi[0] = GetFCNChisquare(pars);
      chi[1] = -nfree;//not working: minimizer->NFree();//need to add nx outside FitKernel
    }
  }

  if(cov){
    minimizer->GetCovMatrix(cov);
  }

  delete minimizer;

 if(fgFitPrintLevel>=0){
    stwa.Stop(); 
    stwa.Print("m"); //"m": milli sec; "u": micro sec
  } 

 //test->
 /*
  for(Int_t ii=0; ii<npar; ii++){
    printf("test0 par %d %e\n", ii, pars[ii]);
    printf("test1 err %d %e\n", ii, errs[ii]);
  }
 */
  //<-test

 return kfail;
}

Int_t XGLUtils::TMinuitFitKernel(const Int_t kml,const Int_t npar, Double_t pars[], Double_t errs[], Double_t chi[], const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs)
{
  //
  //fitkernel for both ML and chi2 fit
  //

  TStopwatch stwa;
  if(fgFitPrintLevel>=0){
    stwa.Start(kTRUE);
  }
  //====>

  TMinuit * mnt = new TMinuit(npar);
  switch(kml){
  case 0:
    mnt->SetFCN(MinFCNChisquare);
    break;
  case 1:
    mnt->SetFCN(MinFCNMaximumLikelihood);
    break;
  case 2:
    mnt->SetFCN(MinFCNBinnedLikelihood);
    break;
  default:
    printf("XGLUtils::TMinuitFitKernel bad kml %d\n", kml);exit(1);
  }
  mnt->SetPrintLevel(fgFitPrintLevel);
    
  //http://wwwasdoc.web.cern.ch/wwwasdoc/WWW/minuit/minmain/chapter2_8.html#SECTION0081000000000000000
  //SET ERRordef<up>
  //Sets the value of UP (default value= 1.), defining parameter errors. Minuit defines parameter errors as the change in parameter value required to change the function value by UP. Normally, for chisquared fits UP=1, and for negative log likelihood, UP=0.5. 
  const Double_t errdef = kml? 0.5 : 1;
  mnt->SetErrorDef(errdef);
    
  Int_t nfree = 0;
  for(Int_t ipar=0;ipar<npar;ipar++){
    const Double_t start=pars[ipar];
    //1e-3*par seems too small
    const Double_t ss= TMath::Abs(pars[ipar])>EPSILON ? TMath::Abs(pars[ipar])*6e-1 : 1e-3;//important if limits are set
    //1e-2; //restore to old values for TRD performance /u/xlu/.task/dedx/drawSummary in /u/xlu/.task/dedx/AllDataSummary/hist/draw_UpdatedStyle_Comb
    //
    const Double_t ll=(lows && lows[ipar]>-1e9)?lows[ipar]:0;
    const Double_t hh=(highs && highs[ipar]>-1e9)?highs[ipar]:0;
    
    mnt->DefineParameter(ipar, Form("p%d",ipar), start, ss, ll,hh);
    nfree++;
    if(pfix && pfix[ipar]){
      mnt->FixParameter(ipar);
      nfree--;
    }
  }

  //http://wwwasdoc.web.cern.ch/wwwasdoc/WWW/minuit/minmain/chapter2_8.html#SECTION0081000000000000000
  //MIGrad[maxcalls] [tolerance]
  //Causes minimization of the function by the method of Migrad, the most efficient and complete single method, recommended for general functions (see also MINImize). The minimization produces as a by-product the error matrix of the parameters, which is usually reliable unless warning messages are produced. The optional argument [maxcalls] specifies the (approximate) maximum number of function calls after which the calculation will be stopped even if it has not yet converged. The optional argument [tolerance] specifies required tolerance on the function value at the minimum. The default tolerance is 0.1, and the minimization will stop when the estimated vertical distance to the minimum (EDM) is less than 0.001*[tolerance]*UP (see [SET ERRordef]SET ERR). 

  Double_t arg[2]={fMaxIter, fTol};
  Int_t migradflag = -999;
  mnt->mnexcm("MIGRAD", arg, 2, migradflag);

  //get par right after MIGRAD before HESSE modifies them, initial pars not changed
  Double_t migp[10];
  Double_t dummy=-999;
  for(Int_t ipar=0;ipar<npar;ipar++){
    mnt->GetParameter(ipar,migp[ipar], dummy);
  }

  Int_t hesseflag = -999;
  mnt->mnexcm("HESSE", arg ,1, hesseflag);

  const Int_t kfail = migradflag+hesseflag;
  if(kfail){
    printf("XGLUtils::FitKernel fitting error !! %d %d keep raw par and set err to -999!!!\n", migradflag, hesseflag);
    if(errs){
      for(Int_t ipar=0;ipar<npar;ipar++){
        errs[ipar]=-999;
      }
    }
    if(chi){
      chi[0]=chi[1]=-999;
    }
  }
  else{
    for(Int_t ipar=0;ipar<npar;ipar++){
      if(errs) {
        mnt->GetParameter(ipar, dummy,errs[ipar]);
      }

      pars[ipar]=migp[ipar];
    }
    if(chi){
      chi[0] = GetFCNChisquare(pars);
      if(nfree!=mnt->GetNumFreePars()){
        printf("XGLUtils::TMinuitFitKernel nfree!=mnt->GetNumFreePars() %d %d\n", nfree, mnt->GetNumFreePars());exit(1);
      }
      chi[1] = -nfree;//need to add nx outside FitKernel
    }
  }

  if(cov){
    mnt->mnemat(cov,npar);
  }

  delete mnt;

  //<===
 
  if(fgFitPrintLevel>=0){
    stwa.Stop(); 
    stwa.Print("m"); //"m": milli sec; "u": micro sec
  } 
 
  //test->
  /*
  for(Int_t ii=0; ii<npar; ii++){
    printf("test0 par %d %e\n", ii, pars[ii]);
    printf("test1 err %d %e\n", ii, errs[ii]);
  }
  */
  //<-test

  return kfail;
}

Int_t XGLUtils::MaximumLikelihoodFit(const Int_t nx, const Double_t xdata[], FFunc ffunc, const Int_t npar, Double_t pars[], Double_t errs[], const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs)
{
  //
  //routine for maximum likelihood fit
  //
  //
  if(nx<=0){
    printf("XGLUtils::MaximumLikelihoodFit error nx<=0 %d\n", nx); exit(1);
  }

  fgFitNStep=0;
  fgFitNData = nx;
  fgFitX = new Double_t[nx];
  memcpy(fgFitX, xdata, sizeof(Double_t)*fgFitNData);
  if(fgFitPrintLevel>=0){
    for(Int_t ix=0; ix<fgFitNData; ix++){
      printf("XGLUtils::MaximumLikelihoodFit %d: %e\n", ix, fgFitX[ix]);
    }
  }
  fgFitFunc = ffunc;

  Double_t tmperr[npar];
  const Int_t kfail = FitKernel(1, npar, pars, tmperr, 0x0, pfix, cov, lows, highs);
  if(errs){
    memcpy(errs, tmperr, sizeof(tmperr));
  }

  delete fgFitX;

  return kfail;
}


void XGLUtils::MinFCNMaximumLikelihood(Int_t &npar, Double_t *gin, Double_t &f, Double_t *pars, Int_t iflag)
{
  f=GetFCNMaximumLikelihood(pars);
}

Double_t XGLUtils::GetFCNMaximumLikelihood(const Double_t pars[])
{
  //
  //function to be minimized in ML fit
  //

  //const Int_t ktest = 1000000;//4165;

  Double_t f = 0;

  /*
  if(fgFitNStep>=ktest){
    for(Int_t ii=0; ii<npar; ii++){
      printf("in MinFCN pars %d %e\n", ii, pars[ii]);
    }
  }
  */

  for(Int_t idata=0; idata<fgFitNData; idata++){
    //if(fgFitNStep>=ktest) printf("test %d %e %e %e\n", idata, fgFitX[idata], pars[0], pars[1]);

    const Double_t lh = fgFitFunc(&(fgFitX[idata]), pars);

    const Double_t eps = 1e-100;
    const Double_t verylarge = -TMath::Log(eps);
    if(lh<eps){
      //means this iteration fails, don't break, otherwise all bad par have the same f (=1e10) which leads to no improvement in the minmization and fail the fit!!
      f += verylarge;
    }
    else{
      f += -TMath::Log(lh);
    }

    //if(fgFitNStep>=ktest) printf("in MinFCN idata %d p0 %e p1 %e data %e lh %e f %e\n", idata, pars[0], pars[1], fgFitX[idata], lh, f);
  }

  //if(fgFitNStep>=ktest) printf("fgFitNStep %d f %e\n", fgFitNStep, f);
  fgFitNStep++;

  return f;
}

Int_t XGLUtils::ChisquareFit(const Int_t nx, const Double_t xdata[], const Int_t nbin, FFunc ffunc, const Int_t npar, Double_t pars[], Double_t errs[], Double_t chi[], const Double_t thres, const Bool_t kxerr, const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs)
{
  return BinnedFit(kFALSE, nx, xdata, nbin, ffunc, npar, pars, errs, chi, thres, kxerr, pfix, cov, lows, highs);
}

Int_t XGLUtils::BinnedLikelihoodFit(const Int_t nx, const Double_t xdata[], const Int_t nbin, FFunc ffunc, const Int_t npar, Double_t pars[], Double_t errs[], const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs)
{
  //keep empty bins
  const Double_t thres = -1;
  return BinnedFit(kTRUE, nx, xdata, nbin, ffunc, npar, pars, errs, 0x0, thres, kFALSE, pfix, cov, lows, highs);
}

Int_t XGLUtils::BinnedFit(const Bool_t kml, const Int_t nx, const Double_t xdata[], const Int_t nbin, FFunc ffunc, const Int_t npar, Double_t pars[], Double_t errs[], Double_t chi[], const Double_t thres, const Bool_t kxerr, const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs)
{
  //
  //chi2fit from unbinned X -> TH1 
  //
  Int_t index[nx];
  TMath::Sort(nx, xdata, index, kFALSE);
  const Double_t xmin = xdata[index[0]];
  const Double_t xmax = xdata[index[nx-1]];
  const Double_t dx = (xmax-xmin)/nbin;

  TH1D hh("chisquarefittmphh","",nbin, xmin-dx, xmax+dx);
  for(Int_t ii=0; ii<nx; ii++){
    hh.Fill(xdata[ii]);
  }

  while(hh.GetBinContent(hh.GetMaximumBin())<10 && hh.GetNbinsX()>10){
    hh.Rebin(2);
  }

  Double_t tmperr[npar];
  Double_t tmpch[10];
  const Int_t kfail = BinnedFit(kml, &hh, ffunc, npar, pars, tmperr, tmpch, thres, kxerr, pfix, cov,lows, highs);
  if(errs){
    memcpy(errs, tmperr, sizeof(tmperr));
  }
  if(chi){
    memcpy(chi, tmpch, sizeof(Double_t)*2);
  }
  return kfail;
}

Int_t XGLUtils::ChisquareFit(const TH1 *hh, FFunc ffunc, const Int_t npar, Double_t pars[], Double_t errs[], Double_t chi[], const Double_t thres, const Bool_t kxerr, const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs)
{
  return BinnedFit(kFALSE, hh, ffunc, npar, pars, errs, chi, thres, kxerr, pfix, cov, lows, highs);
}

Int_t XGLUtils::BinnedLikelihoodFit(const TH1 *hh, FFunc ffunc, const Int_t npar, Double_t pars[], Double_t errs[], const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs)
{
  //keep empty bins
  const Double_t thres = -1;
  return BinnedFit(kTRUE, hh, ffunc, npar, pars, errs, 0x0, thres, kFALSE, pfix, cov, lows, highs);
}

Int_t XGLUtils::BinnedFit(const Bool_t kml, const TH1 *hh, FFunc ffunc, const Int_t npar, Double_t pars[], Double_t errs[], Double_t chi[], const Double_t thres, const Bool_t kxerr, const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs)
{
  //
  //chisquare fit from TH1
  //

  const Int_t nmax = hh->GetNbinsX();
  Double_t xdata[nmax];
  Double_t ydata[nmax];
  Double_t xerr[nmax];
  Double_t yerr[nmax];

  Int_t nx=0;
  for(Int_t ii=hh->GetXaxis()->GetFirst(); ii<=hh->GetXaxis()->GetLast(); ii++){
    const Double_t ey = hh->GetBinError(ii);
    const Double_t yy = hh->GetBinContent(ii);
    //keep all for kml==kTRUE
    //very important for BinnedLikelihood!! only keeping empty bins can the result be good
    if(!kml){
      //skip empty bin
      if(ey<EPSILON){
        if(yy>EPSILON){
          printf("XGLUtils::ChisquareFit error! %d %e %e\n", ii, ey, yy); exit(1);
        }
        continue;
      }

      //skip low statstics bin
      if(yy<thres){
        continue;
      }
    }
    const Double_t xx = hh->GetBinCenter(ii);
    const Double_t ex = hh->GetBinWidth(ii)/TMath::Sqrt(12);

    xdata[nx]=xx;
    ydata[nx]=yy;
    xerr[nx]=ex;
    yerr[nx]=ey;
    nx++;
  }

  return BinnedFit(kml, nx, xdata, ydata, kxerr?xerr:0x0, yerr, ffunc, npar, pars, errs, chi, pfix, cov, lows, highs);
}

Int_t XGLUtils::ChisquareFit(const Int_t nx, const Double_t xdata[], const Double_t ydata[], const Double_t *xerr, const Double_t yerr[], FFunc ffunc, const Int_t npar, Double_t pars[], Double_t errs[], Double_t chi[], const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs)
{
  return BinnedFit(kFALSE, nx, xdata, ydata, xerr, yerr, ffunc, npar, pars, errs, chi, pfix, cov, lows, highs);
}

Int_t XGLUtils::BinnedLikelihoodFit(const Int_t nx, const Double_t xdata[], const Double_t ydata[], FFunc ffunc, const Int_t npar, Double_t pars[], Double_t errs[], const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs)
{
  return BinnedFit(kTRUE, nx, xdata, ydata, 0x0, 0x0, ffunc, npar, pars, errs, 0x0, pfix, cov, lows, highs);
}

Int_t XGLUtils::BinnedFit(const Bool_t kml, const Int_t nx, const Double_t xdata[], const Double_t ydata[], const Double_t *xerr, const Double_t yerr[], FFunc ffunc, const Int_t npar, Double_t pars[], Double_t errs[], Double_t chi[], const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs)
{
  //
  //routine for chisquare fit from arrays X and Y
  //
  //
  if(nx<=0){
    printf("XGLUtils::ChisquareFit error nx<=0 %d\n", nx); exit(1);
  }

  fgFitNStep=0;
  fgFitNData = nx;
  fgFitX = new Double_t[nx];
  fgFitY = new Double_t[nx];
  fgFitEx = new Double_t[nx];
  fgFitEy = new Double_t[nx];
  for(Int_t ix=0; ix<fgFitNData; ix++){
    fgFitX[ix]=xdata[ix];
    fgFitY[ix]=ydata[ix];
    fgFitEx[ix]=xerr?xerr[ix]:0;
    fgFitEy[ix]=yerr[ix];
    if(fgFitPrintLevel>=0){
      printf("XGLUtils::ChisquareFit %d: %e %e %e %e\n", ix, fgFitX[ix], fgFitY[ix], fgFitEx[ix], fgFitEy[ix]);
    }
  }
  fgFitFunc = ffunc;

  const Int_t kfail = FitKernel(kml?2:0, npar, pars, errs, chi, pfix, cov, lows, highs);

  if(chi){
    chi[1] += nx;
    printf("\tXGLUtils::ChisquareFit Chi2/NDOF %.1f/%.0f = %.1f\n\n", chi[0], chi[1], chi[1]? chi[0]/chi[1] : -999);
  }

  delete fgFitX;
  delete fgFitY;
  delete fgFitEx;
  delete fgFitEy;

  return kfail;
}

Double_t XGLUtils::GetFCNChisquare(const Double_t pars[])
{
  //
  //chi2
  //

  //const Int_t ktest = 1000000;//207;//4165;

  Double_t f = 0;

  /*
  if(fgFitNStep>=ktest){
    for(Int_t ii=0; ii<npar; ii++){
      printf("in MinFCN pars %d %e\n", ii, pars[ii]);
    }
  }
  */

  for(Int_t idata=0; idata<fgFitNData; idata++){
    //if(fgFitNStep>=ktest) printf("test %d %e %e %e %e %e %e\n", idata, fgFitX[idata], pars[0], pars[1], pars[2], pars[3], pars[4]);

    const Double_t evaly = fgFitFunc(&(fgFitX[idata]), pars);
    const Double_t dy = fgFitY[idata] - evaly;
    
    Double_t dfdx=0;
    //calculated df/dx if necessary
    if( fgFitEx[idata]!=0){
      const Double_t derix0=fgFitX[idata]-fgFitEx[idata];
      const Double_t derix1=fgFitX[idata]+fgFitEx[idata];
      const Double_t deriy0=fgFitFunc(&derix0, pars);
      const Double_t deriy1=fgFitFunc(&derix1, pars);
      dfdx = (deriy1-deriy0)/(derix1-derix0);
    }
    f += dy*dy/(fgFitEy[idata]*fgFitEy[idata] + dfdx*fgFitEx[idata]*dfdx*fgFitEx[idata]);
    
    //if(fgFitNStep>=ktest) printf("in MinFCN idata %d p0 %e p1 %e data %e f %e dy %e ey %e\n", idata, pars[0], pars[1], fgFitX[idata], f, dy, fgFitEy[idata]);
  }

  //if(fgFitNStep>=ktest) printf("fgFitNStep %d f %e\n", fgFitNStep, f);
  fgFitNStep++;

  return f;
}

void XGLUtils::MinFCNChisquare(Int_t &npar, Double_t *gin, Double_t &f, Double_t *pars, Int_t iflag)
{
  //
  //function to be minimized in chi2 fit
  //

  f = GetFCNChisquare(pars);
}

Double_t XGLUtils::GetFCNBinnedLikelihood(const Double_t pars[])
{
  //
  //binned likelihood
  //

  //const Int_t ktest = 0;//1000000;//207;//4165;

  //the normalization is important for the error estimation, but not the mean!
  //the norm should be the same for all bins, so just get it from first non-0 bin
  Double_t norm = 0;
  for(Int_t idata=0; idata<fgFitNData; idata++){
    if(fgFitEy[idata]>EPSILON){
      if(fgFitY[idata]<EPSILON){
        printf("XGLUtils::GetFCNBinnedLikelihood error y<eps, while ey>eps!! %e %e %e", EPSILON,fgFitEy[idata],fgFitY[idata] );exit(1);
      }
      //y = kN, ey = k sqrt(N), y/ey = sqrt(N)
      const Double_t trueCount = TMath::Power(fgFitY[idata]/fgFitEy[idata],2);
      norm = trueCount/fgFitY[idata];
      break;
    }
  }
  if(norm<EPSILON){
    printf("XGLUtils::GetFCNBinnedLikelihood error norm 0 !! %e\n", norm); exit(1);
  }

  Double_t f = 0;

  for(Int_t idata=0; idata<fgFitNData; idata++){
    const Double_t evaly = fgFitFunc(&(fgFitX[idata]), pars);
    const Double_t prob = TMath::Poisson(fgFitY[idata]*norm, evaly*norm);
    f += -TMath::Log( prob>EPSILON ? prob : EPSILON);
    
    //if(fgFitNStep>=ktest) printf("in MinFCN idata %d p0 %e p1 %e p2 %e x %e y %e y0 %e prob %e f %e\n", idata, pars[0], pars[1], pars[2], fgFitX[idata], fgFitY[idata], evaly, prob, f);
  }

  //if(fgFitNStep>=ktest) printf("fgFitNStep %d f %e\n", fgFitNStep, f);
  //if(fgFitNStep>=ktest) printf("test\n"); exit(1);

  fgFitNStep++;
  return f;
}

void XGLUtils::MinFCNBinnedLikelihood(Int_t &npar, Double_t *gin, Double_t &f, Double_t *pars, Int_t iflag)
{
  //
  //function to be minimized in chi2 fit
  //

  f = GetFCNBinnedLikelihood(pars);
}


TH1D * XGLUtils::GetHfit(const TString hname, FFunc func, const Double_t par[], const Double_t xmin, const Double_t xmax, const Bool_t kbinlogx, const Int_t nbins)
{
  TH1D * h1=new TH1D(hname, "", nbins, xmin, xmax);
  if(kbinlogx){
    BinLogX(h1->GetXaxis());
  }
  for(Int_t ii=1; ii<=h1->GetNbinsX(); ii++){
    const Double_t xx = h1->GetBinCenter(ii);
    const Double_t yy = func(&xx, par);
    h1->SetBinContent(ii, yy);
  }
  return h1;
}

//===================================================================
//                           Histograms
//===================================================================
void XGLUtils::PrintHist(const TString tag, const TH1 * hh)
{
  printf("XGLUtils::PrintHist %s\n", tag.Data());
  Double_t sumre = 0;
  Double_t sumy = 0;
  Int_t nn = 0;
  for(Int_t ii=1; ii<=hh->GetNbinsX(); ii++){
    printf("XGLUtils::PrintHist %s %s ii %d x %e y %e ey %e re %e\n", tag.Data(), hh->GetName(), ii,  hh->GetBinCenter(ii),  hh->GetBinContent(ii), hh->GetBinError(ii), hh->GetBinError(ii)/(hh->GetBinContent(ii)+EPSILON));
    
    if( hh->GetBinContent(ii) >EPSILON){
      sumy += hh->GetBinContent(ii);
      sumre += hh->GetBinError(ii)/(hh->GetBinContent(ii)+EPSILON);
      nn++;
    }
  }
  printf("XGLUtils::PrintHist %s %s meany %e meanre %e\n", tag.Data(), hh->GetName(), sumy/nn, sumre/nn);
}

void XGLUtils::PrintGr(const TString tag, const TGraphAsymmErrors * gr)
{
  printf("XGLUtils::PrintGr %s\n", tag.Data());
  Double_t sumre = 0;
  Double_t sumy = 0;
  Int_t nn = 0;
  for(Int_t ii=0; ii<gr->GetN(); ii++){
    printf("XGLUtils::PrintGr %s %s ii %d x %e exl %e exh %e y %e eyl %e eyh %e rel %e reh %e\n", tag.Data(), gr->GetName(), ii,  gr->GetX()[ii], gr->GetEXlow()[ii], gr->GetEXhigh()[ii], gr->GetY()[ii], gr->GetEYlow()[ii], gr->GetEYhigh()[ii], gr->GetEYlow()[ii]/(gr->GetY()[ii]+EPSILON), gr->GetEYhigh()[ii]/(gr->GetY()[ii]+EPSILON));
    
    if(gr->GetY()[ii]>EPSILON){
      sumy += gr->GetY()[ii];
      sumre += gr->GetEYlow()[ii]/(gr->GetY()[ii]+EPSILON);
      nn++;
    }
  }
  printf("XGLUtils::PrintGr %s %s meany %e meanre %e\n", tag.Data(), gr->GetName(), sumy/nn, sumre/nn);
}

TGraphAsymmErrors * XGLUtils::PublicFilterGraph(const TGraphAsymmErrors * rawgr, const Double_t thres)
{
  TGraphAsymmErrors * outgr = (TGraphAsymmErrors *) rawgr->Clone(Form("%sout", rawgr->GetName()));
  outgr->GetXaxis()->SetLimits(rawgr->GetXaxis()->GetXmin(), rawgr->GetXaxis()->GetXmax());
  Int_t jj = 0;
  for(Int_t ii=0; ii<rawgr->GetN(); ii++){
    const Double_t yy=rawgr->GetY()[ii];
    if(yy<thres){
      continue;
    }

    outgr->GetX()[jj]=rawgr->GetX()[ii];
    outgr->GetEXlow()[jj]=rawgr->GetEXlow()[ii];
    outgr->GetEXhigh()[jj]=rawgr->GetEXhigh()[ii];
    outgr->GetY()[jj]=rawgr->GetY()[ii];
    outgr->GetEYlow()[jj]=rawgr->GetEYlow()[ii];
    outgr->GetEYhigh()[jj]=rawgr->GetEYhigh()[ii];
    jj++;
  }
  outgr->Set(jj);

  return outgr;
}

Bool_t XGLUtils::SimpleScaleToRef(const TGraphAsymmErrors * grStat, const TGraphAsymmErrors * grSys, const TGraphAsymmErrors * grRef, TGraphAsymmErrors *& outStat, TGraphAsymmErrors * & outSys, const Bool_t exitOnError)
{

  TGraphAsymmErrors * tmpStat = (TGraphAsymmErrors *) grStat->Clone(Form("ratio%s",grStat->GetName()));
  TGraphAsymmErrors * tmpSys  = 0x0;
  if(grSys){
    tmpSys = (TGraphAsymmErrors *) grSys->Clone(Form("ratio%s",grSys->GetName()));
  }

  //printf("test XGLUtils::SimpleScaleToRef ======= grstat max %f min %f xtit %s ytit %s\n", grStat->GetMaximum(), grStat->GetMinimum(), grStat->GetXaxis()->GetTitle(), grStat->GetYaxis()->GetTitle());

  const Int_t ngr = grSys? 3 : 2;
  const Int_t refAnchor = 0;
  const Int_t statAnchor = 1;
  const TGraphAsymmErrors * grs[]={grRef, grStat, grSys};

  Int_t id0s[ngr];
  Int_t id1s[ngr];
  const Int_t nptmp = GetCommonLimits(ngr, grs, id0s, id1s);
  if(nptmp<=0){
    if(exitOnError){
      printf("XGLUtils::SimpleScaleToRef nptmp<0 exitOnError %d nptmp %d\n", exitOnError, nptmp); exit(1);
    }
    else{
      return kFALSE;
    }
  }
  
  for(Int_t ip=0; ip<tmpStat->GetN(); ip++){
    Double_t refscal = -999;    

    if(ip>=id0s[statAnchor] && ip<=id1s[statAnchor]){
      refscal = grRef->GetY()[id0s[refAnchor]+(ip-id0s[statAnchor])];
    }

    if(refscal<-990 || fabs(refscal)<EPSILON){
      printf("XGLUtils::SimpleScaleToRef interpolation fail at name %s title %s ip %d npoint %d  refscal %e\n", grStat->GetName(), grStat->GetTitle(), ip, grRef->GetN(),  refscal);
      tmpStat->GetY()[ip]=-999;
      tmpStat->GetEYlow()[ip]=-999;
      tmpStat->GetEYhigh()[ip]=-999;
      if(tmpSys){
        tmpSys->GetY()[ip]=-999;
        tmpSys->GetEYlow()[ip]=-999;
        tmpSys->GetEYhigh()[ip]=-999;
      }
    }
    else{
      tmpStat->GetY()[ip]/=refscal;
      tmpStat->GetEYlow()[ip]/=refscal;
      tmpStat->GetEYhigh()[ip]/=refscal;
      if(tmpSys){
        tmpSys->GetY()[ip]/=refscal;
        tmpSys->GetEYlow()[ip]/=refscal;
        tmpSys->GetEYhigh()[ip]/=refscal;
      }
    }
  } 

  outStat =  PublicFilterGraph(tmpStat, -990);
  if(outStat->GetN()==0){
    printf("XGLUtils::SimpleScaleToRef outStat n=0!!\n");
    for(Int_t it=0; it<grRef->GetN(); it++){
      const Double_t xx = grRef->GetX()[it];
      printf("XGLUtils::SimpleScaleToRef %d %e\n", it, xx);
    }
    exit(1);
  }
  delete tmpStat;

  if(tmpSys){
    outSys =  PublicFilterGraph(tmpSys, -990);
    delete tmpSys;
  }

  return kTRUE;
}

void XGLUtils::ScaleGraphPoint(TGraphAsymmErrors * gr, const Int_t ipoint, const Double_t fac)
{
  const Double_t xx = gr->GetX()[ipoint];
  const Double_t yy = gr->GetY()[ipoint];
  const Double_t el = gr->GetEYlow()[ipoint];
  const Double_t eh = gr->GetEYhigh()[ipoint];
  gr->SetPoint(ipoint,xx, yy*fac);
  gr->SetPointEYlow(ipoint, el*fac);
  gr->SetPointEYhigh(ipoint, eh*fac);
}

void XGLUtils::ScaleGraph(TGraphAsymmErrors * gr, const Double_t fac)
{
  for(Int_t ii=0; ii<gr->GetN(); ii++){
    const Double_t xx = gr->GetX()[ii];
    const Double_t yy = gr->GetY()[ii];
    const Double_t el = gr->GetEYlow()[ii];
    const Double_t eh = gr->GetEYhigh()[ii];
    gr->SetPoint(ii,xx, yy*fac);
    gr->SetPointEYlow(ii, el*fac);
    gr->SetPointEYhigh(ii, eh*fac);
  }
}

void XGLUtils::BinLogX(TAxis *axis) 
{
  //
  // Method for the correct logarithmic binning of histograms
  // copied and modified from AliTPCcalibBase

  const Int_t bins = axis->GetNbins();

  const Double_t from = axis->GetXmin();
  const Double_t to = axis->GetXmax();
  if (from<EPSILON) return;
  Double_t *new_bins = new Double_t[bins + 1];

  new_bins[0] = from;
  const Double_t factor = TMath::Power(to/from, 1./bins);

  for (int i = 1; i <= bins; i++) {
   new_bins[i] = factor * new_bins[i-1];
  }
  axis->Set(bins, new_bins);
  delete [] new_bins;
}


TH1D * XGLUtils::GetMajorCluster(const TH1D * hin, TH1D * & hrecord, Int_t & clsid, Double_t & cut)
{
  TH1D * hout = (TH1D*) hin->Clone(Form("%sout", hin->GetName()));

  hrecord=(TH1D*)hin->Clone(Form("%srec", hin->GetName()));

  clsid = -999;
  cut = -999;

  const Int_t oldnb = XGLUtils::GetNbumps(hrecord);  
  if(oldnb!=1){

  //==============================================

  //if 0,1,2, not in
  while(XGLUtils::GetNbumps(hrecord)>2){
    hrecord->Rebin(2);
  }
    
  const Int_t nbump = XGLUtils::GetNbumps(hrecord);
  if(nbump==0){
    printf("XGLUtils::GetMajorCluster no bump at all!! empty!?\n"); exit(1);
  }
  else if(nbump==1){
    //should still be ok, just because no prominant 2nd cluster show up
    /*
      printf("XGLUtils::GetMajorCluster strange nbump==1 while old!=1 %d\n", oldnb); 
      TCanvas *c1=new TCanvas;
      hrecord->Draw();
      c1->Print("getmajorcluster.eps");
      exit(1);
    */
  }
  else{
    Int_t trigger = 0;
    for(Int_t ii=1; ii<=hrecord->GetNbinsX(); ii++){
      if(hrecord->GetBinContent(ii)>EPSILON){
        trigger = 1;
      }
      if(trigger == 1 && hrecord->GetBinContent(ii)<EPSILON){
        cut = hrecord->GetBinCenter(ii);
        break;
      }
    }
    if(cut<-990){
      printf("XGLUtils::GetMajorCluster cut not found! %d\n", trigger); exit(1);
    }
    const Int_t binc = hout->GetXaxis()->FindBin(cut);
    const Double_t n1 = hout->Integral(1, binc-1);
    const Double_t n2 = hout->Integral(binc+1, hout->GetNbinsX());
    if(fabs(n1-n2)<EPSILON){
      printf("XGLUtils::GetMajorCluster both clusters equal size!! cut %e binc %d n1 %e n2 %e\n", cut, binc, n1, n2); exit(1);
    }
    else if(n1>n2){
      clsid = 1;
      //wipe out n2
      for(Int_t ix=binc+1; ix<=hout->GetNbinsX(); ix++){
        hout->SetBinContent(ix,0); hout->SetBinError(ix, 0);
      }
    }
    else{
      clsid = 2;
      //wipe out n1
      for(Int_t ix=1; ix<=binc-1; ix++){
        hout->SetBinContent(ix,0); hout->SetBinError(ix, 0);
      }
    }
    //not working, hout without rebin has more than 1 bumps
    /*
      const Int_t nbout = XGLUtils::GetNbumps(hout);
      if(nbout!=1){
      printf("XGLUtils::GetMajorCluster nbump of hout !=1 %d\n", nbout); exit(1);
      }
    */
  }
  }

  const Double_t nin  = hin->Integral(1,hin->GetNbinsX()); 
  const Double_t nout = hout->Integral(1, hout->GetNbinsX()); hout->SetEntries(nout);
  printf("XGLUtils::GetMajorCluster cluster fractions %s %e , %e %e\n", hin->GetName(), nout/nin, nin, nout);
  return hout;
}


Int_t XGLUtils::GetNbumps(const TH1D * hh)
{
  TString tag;
  for(Int_t ii=1; ii<=hh->GetNbinsX(); ii++){
    if(hh->GetBinContent(ii)>EPSILON){
      tag+="1";
    }
    else{
      tag+="0";
    }
  }

  while(tag.Contains("11")){
    tag.ReplaceAll("11","1");
  }

  Int_t counter=0;
  for(Int_t ii=0; ii<tag.Length(); ii++){
    if(tag[ii]=='1')
      counter++;
  }

  return counter;
}

TGraphAsymmErrors * XGLUtils::MeasurementConversion(const TString hn, const TGraphAsymmErrors * gv1, const TH2D * hy2x1, const Bool_t kcutEdge, const Double_t thres)
{
  //hy2x1 X matches gv1 X
  //hy2x1 Y matches output X

  //example
  /*
//should also work for restrict x range!


TGraphAsymmErrors *g1=new TGraphAsymmErrors(3)
g1->SetPoint(0,4,7)
g1->SetPoint(1,5,3)
g1->SetPoint(2, 6, 2)
g1->SetPointError(0,0,0,1,2)
g1->SetPointError(1,0,0,1,2)
g1->SetPointError(2,0,0,1,2)
g1->Draw("ap*")

TH2D *hy2x1=new TH2D("hh","",7,1.5,8.5,4,99,103)
hy2x1->SetBinContent(1,1,100)
hy2x1->SetBinContent(1,2,100)
hy2x1->SetBinContent(1,3,100)
hy2x1->SetBinContent(1,4,100)
hy2x1->SetBinContent(2,1,100)
hy2x1->SetBinContent(2,2,100)
hy2x1->SetBinContent(2,3,100)
hy2x1->SetBinContent(2,4,100)
hy2x1->SetBinContent(3,1,3)
hy2x1->SetBinContent(3,2,9)
hy2x1->SetBinContent(3,3,2)
hy2x1->SetBinContent(3,4,0)
hy2x1->SetBinContent(4,1,8)
hy2x1->SetBinContent(4,2,5)
hy2x1->SetBinContent(4,3,6)
hy2x1->SetBinContent(4,4,1)
hy2x1->SetBinContent(5,1,1)
hy2x1->SetBinContent(5,2,5)
hy2x1->SetBinContent(5,3,4)
hy2x1->SetBinContent(5,4,7)
hy2x1->SetBinContent(6,1,100)
hy2x1->SetBinContent(6,2,100)
hy2x1->SetBinContent(6,3,100)
hy2x1->SetBinContent(6,4,100)
hy2x1->SetBinContent(7,1,100)
hy2x1->SetBinContent(7,2,100)
hy2x1->SetBinContent(7,3,100)
hy2x1->SetBinContent(7,4,100)
hy2x1->Draw("text")

hy2x1->GetXaxis()->SetRangeUser(g1->GetX()[0], g1->GetX()[g1->GetN()-1]);


gout=XGLUtils::MeasurementConversion("test",g1, hy2x1)
for(int ii=0; ii<gout->GetN(); ii++){printf("%d/%d %e %e %e -- %e %e %e\n", ii, gout->GetN(), gout->GetX()[ii], gout->GetEXlow()[ii], gout->GetEXhigh()[ii], gout->GetY()[ii], gout->GetEYlow()[ii], gout->GetEYhigh()[ii]);}

//===

should expect gout
X: 99.5, 100.5, 101.5, 102.5, all errors 0.5
Y: 3.917, 4.632, 3.333, 2.125,
EYlow: 0.717, 0.602, 0.624, 0.884
EYhigh: 2*EYlow

   */

  const Int_t matrixFirstX = hy2x1->GetXaxis()->GetFirst();
  //check nbins
  const Int_t nv1 = gv1->GetN();
  const Int_t matrixNx = hy2x1->GetXaxis()->GetLast() - matrixFirstX + 1;
  if(matrixNx!=nv1){
    printf("XGLUtils::MeasurementConversion hy2x1 nv1 != hy2xinx %s: %d , %s %d\n", gv1->GetName(), nv1, hy2x1->GetName(), matrixNx); return 0x0;
  }

  //check v1 binning
  for(Int_t ii=0; ii<nv1; ii++){
    if( fabs( gv1->GetX()[ii] - hy2x1->GetXaxis()->GetBinCenter(ii+matrixFirstX) ) > EPSILON ){
      printf("XGLUtils::MeasurementConversion gv1 hy2x1 binning does NOT match at ii %d gv1 %e hy2x1 %e\n", ii, gv1->GetX()[ii], hy2x1->GetXaxis()->GetBinCenter(ii+matrixFirstX)); return 0x0;
    }
  }

  const Int_t nv2 = hy2x1->GetNbinsY();
  TGraphAsymmErrors * gout = new TGraphAsymmErrors;
  gout->SetName(hn);
  gout->SetTitle(Form("converted from %s %s", gv1->GetName(), hy2x1->GetName()));

  for(Int_t i2 = 0; i2< nv2; i2++){
    Double_t ff = 0;
    Double_t sl2 = 0;
    Double_t sh2 = 0;

    //Double_t Integral(Int_t binx1, Int_t binx2, Int_t biny1, Int_t biny2, Option_t* option = "") const 
    const Double_t nor = hy2x1->Integral(matrixFirstX, nv1+matrixFirstX-1, i2+1, i2+1);

    //if on the edge there is counts, it means the acceptance is not complete
    const Double_t edgeThres = nor*1e-2;

    //the bin is large enough, no need to have two bins
    //const Double_t edgeCount = hy2x1->Integral(matrixFirstX, matrixFirstX+1, i2+1, i2+1) + hy2x1->Integral(nv1+matrixFirstX-2, nv1+matrixFirstX-1, i2+1, i2+1);

    const Double_t edgeCount = hy2x1->GetBinContent(matrixFirstX,i2+1) + hy2x1->GetBinContent(nv1+matrixFirstX-1, i2+1);

    if(nor>thres){
      if(!kcutEdge || edgeCount< edgeThres){

        for(Int_t i1 = 0; i1 < nv1; i1++){
          const Double_t v1f = gv1->GetY()[i1];
          const Double_t ww = hy2x1->GetBinContent(i1+matrixFirstX, i2+1)/nor;

          if(v1f<EPSILON || ww<EPSILON){
            continue;
          }

          const Double_t v1el = gv1->GetEYlow()[i1];
          const Double_t v1eh = gv1->GetEYhigh()[i1];
          
          const Double_t dw = TMath::Sqrt(ww*(1-ww)/nor);

          //v1elh correspond to the whole pbin ~1/sqrt(totPstat)
          //(p,z) corresponds to only bin content ~1/sqrt(bincontent)
          const Double_t totPstat = hy2x1->Integral(i1+matrixFirstX, i1+matrixFirstX, 0, nv2);
          const Double_t inflat = totPstat/hy2x1->GetBinContent(i1+matrixFirstX, i2+1);

          ff += v1f*ww;
          //assuming uncorrelated errors
          sl2 += TMath::Power(v1f*ww, 2) * ( TMath::Power(dw/ww,2) + TMath::Power(v1el/v1f,2)*inflat );//ww*ww*v1el*v1el; ww*ww*v1eh*v1eh;
          sh2 += TMath::Power(v1f*ww, 2) * ( TMath::Power(dw/ww,2) + TMath::Power(v1eh/v1f,2)*inflat );//ww*ww*v1el*v1el; ww*ww*v1eh*v1eh;
          
          //printf("XGLUtils::MeasurementConversion test %s i2 %d/%d yy %e nor %e , i1 %d/%d, v1f %e v1el %e v1eh %e ww*nor %e ff %e sl2 %e sh2 %e\n", gout->GetName(), i2, nv2, hy2x1->GetYaxis()->GetBinCenter(i2+1),  nor, i1, nv1, v1f, v1el, v1eh, ww*nor, ff, sl2, sh2);
        }

        const Double_t efl = TMath::Sqrt(sl2);
        const Double_t efh = TMath::Sqrt(sh2);

        const Int_t np = gout->GetN();          
        gout->SetPoint(np, hy2x1->GetYaxis()->GetBinCenter(i2+1), ff);
        gout->SetPointError(np, hy2x1->GetYaxis()->GetBinWidth(i2+1)/2, hy2x1->GetYaxis()->GetBinWidth(i2+1)/2, efl, efh);

      }
      else{
        printf("XGLUtils::MeasurementConversion not calculated %s i2 %d nor %e kcutEdge %d edgeCount %e %e %e edgeThres %e\n", gout->GetName(), i2, nor, kcutEdge, hy2x1->GetBinContent(matrixFirstX,i2+1), hy2x1->GetBinContent(nv1+matrixFirstX-1, i2+1), edgeCount, edgeThres);
        //do not record any more if reaching edge
        if(gout->GetN()){
          break;
        }
      }
    }


  }

  return gout;
}

TGraphAsymmErrors * XGLUtils::GetPIDCont(TH1D *h1, TH1D *h2)
{
  TH1D * hhs[]={h1, h2};

  TGraphAsymmErrors * pid = new TGraphAsymmErrors;
  const Int_t lastbin = hhs[0]->GetNbinsX()+1;
  for(Int_t ibin=1; ibin < lastbin; ibin++){
    //eff for electron
    const Double_t eff_nu = hhs[0]->Integral(ibin, lastbin);
    const Double_t eff_de = hhs[0]->Integral(0, lastbin);
    const Double_t eff_va = eff_nu/eff_de;
    const Double_t eff_er = TMath::Sqrt( eff_va*(1-eff_va)/(eff_de+EPSILON) + EPSILON );


    //absolute contamination for pion in electron
    const Double_t cont_nu = hhs[1]->Integral(ibin, lastbin);
    const Double_t cont_de = eff_nu;
    const Double_t cont_va = cont_nu/cont_de;
    const Double_t cont_er = cont_va * TMath::Sqrt( 1/(cont_nu+EPSILON) + 1/(cont_de+EPSILON) );

    const Int_t ip = pid->GetN();
    pid->SetPoint(ip, eff_va, cont_va);
    pid->SetPointError(ip, eff_er, eff_er, cont_er, cont_er);
  }
  return pid;

}

TList * XGLUtils::Get2DPIDEff(TH2D *h1, TH2D *h2, const Int_t nxcut, const Double_t xcuts[])
{
  TH2D *hhs[]={h1,h2};

  TList * ll = new TList;

  const Int_t lastx = hhs[0]->GetNbinsX()+1;
  for(Int_t icut=0; icut<nxcut; icut++){
    const Int_t ix = hhs[0]->GetXaxis()->FindBin(xcuts[icut]);
    const Double_t realcut = hhs[0]->GetXaxis()->GetBinLowEdge(ix);

    TGraphAsymmErrors * igr = new TGraphAsymmErrors;
    igr->SetName(Form("gr%d", icut));
    igr->SetTitle(Form("cut x>%e", realcut));
    ll->Add(igr);

    const Int_t lasty = hhs[0]->GetNbinsY()+1;
    for(Int_t iy= 1; iy< lasty; iy++){

      Double_t effs[]={-999,-999};
      Double_t errs[]={-999,-999};
      for(Int_t ih=0; ih<2; ih++){
   
        const Double_t nu  = hhs[ih]->Integral(ix, lastx, iy, lasty);
        const Double_t den = hhs[ih]->Integral(0, lastx, 0, lasty);
        const Double_t eff = nu/den;
        const Double_t err = TMath::Sqrt( eff*(1-eff)/(den+EPSILON) + EPSILON );
        effs[ih]=eff;
        errs[ih]=err;

      }

      const Int_t ip = igr->GetN();
      igr->SetPoint(ip, effs[0], effs[1]);
      igr->SetPointError(ip, errs[0],errs[0],errs[1],errs[1]);
    }
  }

  printf("\n\n===================================================\nTList * XGLUtils::Get2DPIDEff list:\n\n");
  ll->ls();
  printf("=======================================================\n\n");
  return ll;
}

TGraphAsymmErrors * XGLUtils::GetPIDEff(TH1D *h1, TH1D *h2)
{
  TH1D * hhs[]={h1, h2};
  
  /*
  const Double_t na = ha->GetNbinsX();
  const Double_t xamin = ha->GetXaxis()->GetBinLowEdge(1);
  const Double_t xamax = ha->GetXaxis()->GetBinUpEdge(na);

  const Double_t nb = hb->GetNbinsX();
  const Double_t xbmin = hb->GetXaxis()->GetBinLowEdge(1);
  const Double_t xbmax = hb->GetXaxis()->GetBinUpEdge(nb);

  if(na!=nb || fabs(xamin-xbmin)>EPSILON || fabs(xamax-xbmax)>EPSILON){
    printf("TGraphAsymmErrors * XGLUtils::GetPIDEff ha hb do not match! %f %f %f %f %f %f\n", na, nb, xamin, xbmin, xamax, xbmax); exit(1);
  }
  */

  TGraphAsymmErrors * pid = new TGraphAsymmErrors;
  const Int_t lastbin = hhs[0]->GetNbinsX()+1;
  for(Int_t ibin=1; ibin < lastbin; ibin++){
    Double_t effs[]={-999,-999};
    Double_t errs[]={-999,-999};
    for(Int_t ih=0; ih<2; ih++){
      const Double_t nu  = hhs[ih]->Integral(ibin, lastbin);
      const Double_t den = hhs[ih]->Integral(0, lastbin);
      const Double_t eff = nu/den;
      const Double_t err = TMath::Sqrt( eff*(1-eff)/(den+EPSILON) + EPSILON );
      effs[ih]=eff;
      errs[ih]=err;
    }
    const Int_t ip = pid->GetN();
    pid->SetPoint(ip, effs[0], effs[1]);
    pid->SetPointError(ip, errs[0],errs[0],errs[1],errs[1]);
  }
  return pid;
}

void XGLUtils::GetGrSysMode1(TList *ll, const TString tag)
{
  printf("XGLUtils::GetGrSysMode1 over-estimate if all components are correlated. Use mode0 instead!!\n"); exit(1);

  TGraphAsymmErrors * raw0 = (TGraphAsymmErrors*) ll->At(0);
  if(!raw0){
    printf("XGLUtils::GetSysGraph no raw0!!\n"); ll->ls(); exit(1);
  }

  TGraphAsymmErrors * grstat = (TGraphAsymmErrors*) raw0->Clone(Form("grstat%s", tag.Data()));
  TGraphAsymmErrors * grsys  = (TGraphAsymmErrors*) raw0->Clone(Form("grsys%s", tag.Data()));

  grsys->SetFillStyle(3001);
  grsys->SetFillColor(kBlack);

  //==================================================
  //must get ngraph before adding any new graph!!
  ll->ls();
  const Int_t ng = ll->GetEntries();
  ll->Add(grstat);
  ll->Add(grsys);

  const Int_t npoint = grsys->GetN();
  printf("XGLUtils::GetSysGraph ng %d np %d\n", ng, npoint);

  Double_t xx[npoint];
  Double_t yy[npoint];
  Double_t sumw2[npoint];

  for(Int_t ii=0; ii<npoint; ii++){
    Double_t tmpx=-999, tmpy=-999;
    grsys->GetPoint(ii, tmpx, tmpy);
    xx[ii]=tmpx;
    yy[ii]=tmpy;

    printf("XGLUtils::GetSysGraph reading %s gr0 %d %e %e\n", tag.Data(), ii, tmpx, tmpy);

    sumw2[ii]=0.0;
  }

  TGraphAsymmErrors * gdeltas[ng];
  
  //start from gr1
  for(Int_t ig=1; ig<ng; ig++){
    TGraphAsymmErrors * gr = (TGraphAsymmErrors*) ll->At(ig);
    if(!gr){
      printf("XGLUtils::GetSysGraph no gr %d!!\n", ig); ll->ls(); exit(1);
    }

    if(gr->GetN()!=npoint){
      printf("XGLUtils::GetSysGraph gr->GetN()!=npoint %d %d\n", gr->GetN(), npoint); exit(1);
    }

    gdeltas[ig] = new TGraphAsymmErrors;
    gdeltas[ig]->SetName(Form("gdelta%sf%d", tag.Data(),ig));
    ll->Add(gdeltas[ig]);

    for(Int_t ix=0; ix<npoint; ix++){
      Double_t tmpx=-999, tmpy=-999;
      gr->GetPoint(ix, tmpx, tmpy);
      printf("XGLUtils::GetSysGraph reading %s %d %e %e\n", gr->GetName(), ix, tmpx, tmpy);

      const Double_t eps = EPSILON;
      if(TMath::Abs(xx[ix]-tmpx)>eps){
        printf("XGLUtils::GetSysGraph xx[ix]!=tmpx %e %e %e %e\n", xx[ix], tmpx, TMath::Abs(xx[ix]-tmpx), eps);exit(1);
      }

      const Double_t dd = tmpy-yy[ix];
      gdeltas[ig]->SetPoint(ix, tmpx, dd);
      sumw2[ix] += dd*dd;
      printf("XGLUtils::GetSysGraph ig %d ix %d dd %e sumw2 %e\n", ig, ix, dd, sumw2[ix]);
    }
  }

  TGraphAsymmErrors * grcontrs[ng];
  for(Int_t ig=0; ig<ng; ig++){
    grcontrs[ig]=0x0;
  }

  for(Int_t ii=0; ii<npoint; ii++){
    const Double_t ey = TMath::Sqrt(sumw2[ii]);
    printf("XGLUtils::GetSysGraph ii %d sumw2 %e ey %e\n",ii, sumw2[ii], ey);

    for(Int_t ig=1; ig<ng; ig++){
      Double_t dy=-999, tmpx=-999;
      gdeltas[ig]->GetPoint(ii,tmpx, dy);
      const Double_t contr = TMath::Power(dy/ey,2);
      printf("XGLUtils::GetSysGraph tag %s ii %d ig %d ey %e dy %e contribution %e\n", tag.Data(), ii, ig, ey, dy, contr);
      
      if(!grcontrs[ig]){
        grcontrs[ig] = new TGraphAsymmErrors;
        grcontrs[ig]->SetName(Form("grcontr%sf%d", tag.Data(), ig));
        ll->Add(grcontrs[ig]);
      }

      grcontrs[ig]->SetPoint(ii, tmpx, contr);
    }

    grsys->SetPointEYlow(ii,ey);
    grsys->SetPointEYhigh(ii,ey);
  }
}

/*
void XGLUtils::GetGrChi2(TList *ll, Double_t grchi2[])
{
  //Chi2 here actually is the deviation: |y/mean-1|

  const Int_t ng = ll->GetEntries();
  TGraphAsymmErrors * grs[ng];
  for(Int_t igr=0; igr<ng; igr++){
    grs[igr] = (TGraphAsymmErrors*) ll->At(igr);
    if(!grs[igr]){
      printf("XGLUtils::GetGrChi2 no gr %d!!\n", igr); ll->ls(); exit(1);
    }
  }

  const Int_t npoint = grs[0]->GetN();
  printf("XGLUtils::GetGrChi2 ng %d np %d\n", ng, npoint);
  ll->ls();

  //clean up first by cutting on the chi2 w.r.t. the mean of all graph
  //assuming all graphs each has npoint
  //assuming the x's are all the same in different graphs, already checked later explicitly
  Double_t meany[npoint];
  for(Int_t ipo=0; ipo<npoint; ipo++){
    meany[ipo]=0;

    Double_t mw = 0;
    for(Int_t igr=0; igr<ng; igr++){
      //const Double_t err = 0.5*(grs[igr]->GetEYlow()[ipo]+grs[igr]->GetEYhigh()[ipo]);
      const Double_t iw = 1; //1/(err*err+EPSILON);
      meany[ipo] += grs[igr]->GetY()[ipo] * iw;
      mw+=iw;
    }

    meany[ipo] /= mw;
  }

  Double_t chi2[ng];
  //no need to get nsigma: suppose all point have the same but large sigma, then it is still 0-sigma chi2, but it should in fact be rejected!
  //Double_t meanGraphChi2 = 0;
  //Double_t sum2 = 0;
  for(Int_t igr=0; igr<ng; igr++){
    chi2[igr] = 0;

    Int_t icount=0;
    for(Int_t ipo=0; ipo<npoint; ipo++){
      if(meany[ipo]>EPSILON){
        //const Double_t err = 0.5*(grs[igr]->GetEYlow()[ipo]+grs[igr]>GetEYhigh()[ipo]);

        //simply the deviation
        //no need to divide mean if later all particles are summed
        //so now simply the difference
        const Double_t deno = 1;//meany[ipo];//(err+EPSILON);
        chi2[igr] += TMath::Abs( (grs[igr]->GetY()[ipo]-meany[ipo])/deno ); //TMath::Power( (grs[igr]>GetY()[ipo]-meany[ipo])/deno, 2);
        icount++;
      }
    }
    if(icount==0){
      printf("XGLUtils::GetGrChi2 bad chi2 icount!! %d\n", icount); exit(1);
    }
    chi2[igr]/=icount;

    //meanGraphChi2 += chi2[igr];
    //sum2 += TMath::Power(chi2[igr],2);
  }
  //meanGraphChi2 /= ng;
  //not good to use, because sigma can be 0 if ng=1 then nsigma = inf
  //const Double_t sigmaGraphChi2 = TMath::Sqrt(sum2/ng-meanGraphChi2*meanGraphChi2);
  
  for(Int_t igr=0; igr<ng; igr++){
    //simply the deviation
    //const Double_t nsigchi2 = TMath::Abs(chi2[igr]-meanGraphChi2)/meanGraphChi2;   //sigmaGraphChi2;
    grchi2[igr]=chi2[igr]; //nsigchi2;
    //printf("XGLUtils::GetGrChi2 chi2 gr %d: %e %e - %e %e -- %d %d\n", igr, chi2[igr], nsigchi2, meanGraphChi2, sigmaGraphChi2, ng, npoint);
    printf("XGLUtils::GetGrChi2 ng %d, npoint %d, chi2[%d] %e\n", ng, npoint, igr, chi2[igr]);
  }
}
*/


Bool_t XGLUtils::GetGrPointLimits(const TGraphAsymmErrors * gr, Double_t & gmin, Double_t & gmax, Bool_t & klogx)
{
  if(!gr){
    printf("XGLUtils::GetGrPointLimits gr null\n"); exit(1);
  }
  if(gr->GetN()==0){
    printf("XGLUtils::GetGrPointLimits gr n=0!!\n"); exit(1);
  }

  gmin = 1e10;
  gmax = -1e10;
  for(Int_t ii=0; ii<gr->GetN(); ii++){
    const Double_t xx = gr->GetX()[ii];
    if(xx>gmax)
      gmax = xx;
    if(xx<gmin)
      gmin = xx;
  }

  const Double_t e0 = gr->GetEXlow()[0];
  const Double_t e1 = gr->GetEXhigh()[0];

  if(fabs(e0-e1)>EPSILON)
    klogx = kTRUE;
  else
    klogx = kFALSE;

  Bool_t kup=kTRUE;
  if( fabs(gr->GetX()[0]-gmin)>EPSILON && fabs(gr->GetX()[gr->GetN()-1]-gmax)<EPSILON){
    kup = kTRUE;
  }
  else{
    kup = kFALSE;
  }

  return kup;
}


Int_t XGLUtils::GetCommonLimits(const Int_t ngr,  const TGraphAsymmErrors * grs[], Int_t id0s[], Int_t id1s[])
{
  Double_t maxlowerpoint= -1e2;
  Double_t minupperpoint = 1e2;
  for(Int_t ii=0; ii<ngr; ii++){
    Double_t gmin = -999;
    Double_t gmax = -999;
    Bool_t klogx = kFALSE;
    if(!grs[ii]->GetN()){
      printf("XGLUtils::GetCommonLimits::bad graph n=0 igr %d %s %s\n", ii, grs[ii]->GetName(), grs[ii]->GetTitle()); 
      return 0;
      //exit(1);
    }
    GetGrPointLimits(grs[ii], gmin, gmax, klogx);
    if(gmin>maxlowerpoint){
      maxlowerpoint = gmin;
    }
    if(gmax<minupperpoint){
      minupperpoint = gmax;
    }
    printf("XGLUtils::GetCommonLimits::loop in graph igr %d/%d min %e max %e maxlowerpoint %e minupperpoint %e\n", ii, ngr, gmin, gmax, maxlowerpoint, minupperpoint);
  }

  Int_t npoint = -999;
  for(Int_t igr=0; igr<ngr; igr++){
    id0s[igr]=-999;
    id1s[igr]=-999;
    for(Int_t pp=0; pp<grs[igr]->GetN(); pp++){
      const Double_t xx = grs[igr]->GetX()[pp];
      if(fabs(xx-maxlowerpoint)<EPSILON){
        id0s[igr] = pp;
      }
      if(fabs(xx-minupperpoint)<EPSILON){
        id1s[igr] = pp;
      }
    }
    if(id0s[igr]<0){
      printf("XGLUtils::GetCommonLimits::not found for %s igr %d maxlowerpoint %e , %e %e\n", grs[igr]->GetName(), igr, maxlowerpoint, grs[igr]->GetX()[0], grs[igr]->GetX()[0]-maxlowerpoint); exit(1);
    }
    if(id1s[igr]<0){
      printf("XGLUtils::GetCommonLimits::not found for igr %d minupperpoint %e\n", igr, minupperpoint); 
      for(Int_t pp=0; pp<grs[igr]->GetN(); pp++){
        const Double_t xx = grs[igr]->GetX()[pp];
        printf("pp %d xx %e del %e\n", pp, xx, xx-minupperpoint);
      }
      exit(1);
    }
    printf("XGLUtils::GetCommonLimits::results igr %d start %d maxlowerpoint %e end %d minupperpoint %e\n", igr, id0s[igr], maxlowerpoint, id1s[igr], minupperpoint);

    const Int_t tmpnp = id1s[igr]-id0s[igr]+1;
    if(npoint<0){
      npoint = tmpnp;
    }
    else if(tmpnp!=npoint){
      printf("XGLUtils::GetCommonLimits::npoint not stable!! igr %d %d %d\n", igr, tmpnp, npoint); 
      return -999;
      //exit(1);
    }
  }
  return npoint;
}


void XGLUtils::GetGrSysMode0(TList *ll, const TString tag, const Bool_t *badgr, const Bool_t kfit, const Double_t hmin, const Double_t hmax)
{
  printf("\n\n********************** XGLUtils::GetGrSysMode0 ********************** \n\n");
  ll->ls();

  const Int_t rawn = ll->GetEntries();

  const TGraphAsymmErrors * grall[rawn];
  Int_t id0s[rawn];
  Int_t id1s[rawn];
  for(Int_t ii=0; ii<rawn; ii++){
    grall[ii] = 0x0;
    id0s[ii]=id1s[ii]=-999;
  }
  Int_t ng = 0;
  for(Int_t ii=0; ii<rawn; ii++){
    if(badgr && badgr[ii]){
      //skip bad gr
      printf("XGLUtils::GetGrSysMode0 %s gr %d finally-kicked-out!\n", tag.Data(), ii);
      continue;
    }
    grall[ng]=(TGraphAsymmErrors*) ll->At(ii);
    ng++;
  }

  const Int_t npoint = GetCommonLimits(ng, grall, id0s, id1s);
  if(npoint<0){
    printf("XGLUtils::GetGrSysMode0 npoint<0 %d GetCommonLimits not stable?\n", npoint); exit(1);
  }
  printf("XGLUtils::GetGrSysMode0 ng %d np %d\n", ng, npoint);

  //output to ll by adding in the very end, not to disturb ng

  TGraphAsymmErrors * grstat = new TGraphAsymmErrors;
  grstat->SetName(Form("grstat%s", tag.Data()));

  TGraphAsymmErrors * grsys = new TGraphAsymmErrors;
  grsys->SetName(Form("grsys%s", tag.Data()));
  grsys->SetFillStyle(3004);
  grsys->SetFillColor(kGray);

  ll->Add(grstat);
  ll->Add(grsys);

  if(!npoint){
    return;
  }

  //==================================================

  Double_t xx[npoint];
  Double_t exl[npoint];
  Double_t exh[npoint];
  Double_t yy[npoint];
  Double_t sumw2y[npoint];

  Double_t eyl[npoint];
  Double_t eyh[npoint];
  Double_t ymin[npoint];
  Double_t ymax[npoint];
  Int_t gmin[npoint];
  Int_t gmax[npoint];
  TH1D *hdists[npoint];

  TH2D * h2d=new TH2D(Form("h2d%s",tag.Data()),"",50, 0.1, 100, 100, hmin, hmax); ll->Add(h2d);
  BinLogX(h2d->GetXaxis());

  Int_t nread = 0;

  const Int_t hnbx = 2000;

  for(Int_t ix=0; ix<npoint; ix++){
    xx[ix]=-999;
    exl[ix]=0;
    exh[ix]=0;
    yy[ix]=0;
    sumw2y[ix]=0;
    eyl[ix]=0;
    eyh[ix]=0;

    ymin[ix]= 1e10;
    ymax[ix]=-1e10;
    gmin[ix]=-1;
    gmax[ix]=-1;

    hdists[ix] = 0x0;

    //=============================
    //loop over graphs
    Double_t grstaterrl[ng];
    Double_t grstaterrh[ng];

    for(Int_t ig=0; ig<ng; ig++){
      grstaterrl[ig] = 1E10;
      grstaterrh[ig] = 1E10;

      const TGraphAsymmErrors * gr = grall[ig];

      if(gr->GetN()<npoint){
        printf("XGLUtils::GetGrSysMode0 ig %d/%d gr->GetN()<npoint %d %d\n", ig, ng, gr->GetN(), npoint); exit(1);
      }

      const Int_t iraw = id0s[ig]+ix;

      Double_t tmpx=-999, tmpy=-999;
      gr->GetPoint(iraw, tmpx, tmpy);
      printf("XGLUtils::GetGrSysMode0 reading %s %d %e %e\n", gr->GetName(), iraw, tmpx, tmpy);

      if(xx[ix]<-990){
        xx[ix]=tmpx;
        exl[ix]=gr->GetEXlow()[iraw];
        exh[ix]=gr->GetEXhigh()[iraw];
      }
      else{
        const Double_t eps = EPSILON;
        if(TMath::Abs(xx[ix]-tmpx)>eps){
          printf("XGLUtils::GetGrSysMode0 xx[%d]!=tmpx %d, %e %e %e %e\n", ix, iraw, xx[ix], tmpx, TMath::Abs(xx[ix]-tmpx), eps);exit(1);
        }
      }

      if(!hdists[ix]){
        hdists[ix] = new TH1D(Form("hdist%s_%d", tag.Data(), ix), Form("xx %e", tmpx), hnbx, hmin, hmax);    ll->Add(hdists[ix]);
      }

      h2d->Fill(tmpx, tmpy);
      hdists[ix]->Fill(tmpy);

      yy[ix]+=tmpy;
      sumw2y[ix]+=tmpy*tmpy;

      //eyl[ix]+= gr->GetEYlow()[iraw];
      //eyh[ix]+= gr->GetEYhigh()[iraw];

      grstaterrl[ig] = gr->GetEYlow()[iraw];
      grstaterrh[ig] = gr->GetEYhigh()[iraw];

      if(tmpy>ymax[ix]){
        ymax[ix]=tmpy;
        gmax[ix]=ig;
      }

      if(tmpy<ymin[ix]){
        ymin[ix]=tmpy;
        gmin[ix]=ig;
      }

      if(ix==0){
        nread++;
      }
    }

    //Double_t XGLUtils::TruncatedMean(const Int_t nx, const Double_t xdata[], const Double_t lowfrac, const Double_t highfrac, Double_t * grms, Double_t * gerr, Double_t *wws)
    //printf("testtest ng %d \n", ng);
    for(Int_t itest=0; itest<ng; itest++){
      if(TMath::IsNaN(grstaterrl[itest])) grstaterrl[itest]=1e3;
      if(TMath::IsNaN(grstaterrh[itest])) grstaterrh[itest]=1e3;

      //printf("testtest %d %e %e\n", itest, grstaterrl[itest], grstaterrh[itest]);
    }
    const Double_t lowf = 0.05;
    const Double_t highf = 0.2;

    const Int_t itmpstart = (Int_t)(lowf*ng);
    const Int_t itmpstop  = (Int_t)(highf*ng);
    if(itmpstop-itmpstart>0){
      eyl[ix] = TruncatedMean(ng, grstaterrl, lowf, highf);
      eyh[ix] = TruncatedMean(ng, grstaterrh, lowf, highf);
    }
    else{
      if(ng!=1){
        printf("XGLUtils::GetGrSysMode0 wrong ng lowf highf %d %e %e\n", ng, lowf, highf);exit(1);
      }
      else{
        eyl[ix]= grstaterrl[0];
        eyh[ix]= grstaterrh[0];
      }
    }
  }

  for(Int_t ii=0; ii<npoint; ii++){

    Double_t pars[10]={ hdists[ii]->Integral(0,hdists[ii]->GetNbinsX()+1)*hdists[ii]->GetBinWidth(1), hdists[ii]->GetMean(), hdists[ii]->GetRMS() };

    if(kfit){
      /*
      printf("XGLUtils::GetGrSysMode0 should not be using kfit!!\n"); exit(1);

      const Double_t minpeak = 20;
      while( hdists[ii]->GetBinContent(hdists[ii]->GetMaximumBin())< minpeak ){
        hdists[ii]->Rebin(2);
        if(hdists[ii]->GetNbinsX()<5){
          printf("XGLUtils::GetGrSysMode0 no peak found! %d %d %e %d\n", hdists[ii]->GetNbinsX(), hdists[ii]->GetMaximumBin(), hdists[ii]->GetBinContent(hdists[ii]->GetMaximumBin()), nread); exit(1);
        }
      }
      
      FFunc ff=RawGaus;
      const Int_t npar = 3;
      Double_t errs[10], chi[10];

      TH1D *tmph=(TH1D*) hdists[ii]->Clone(Form("%stmp", hdists[ii]->GetName())); ll->Add(tmph);
      const Int_t nsig = 3;
      const Double_t tmpxmin = tmph->GetMean()-nsig*tmph->GetRMS();
      const Double_t tmpxmax = tmph->GetMean()+nsig*tmph->GetRMS();
      for(Int_t itmp=0; itmp<=tmph->GetNbinsX()+1; itmp++){
        const Double_t tmpxx = tmph->GetBinCenter(itmp);
        if(tmpxx<tmpxmin || tmpxx>tmpxmax){
          tmph->SetBinContent(itmp,0);
        }
      }

      pars[1]=tmph->GetMean();
      pars[2]=tmph->GetRMS();

      const Int_t kfail = ChisquareFit(tmph, ff, npar, pars, errs, chi);

      TH1D *hfit = GetHfit(Form("%sfit", hdists[ii]->GetName()), ff, pars, hdists[ii]->GetXaxis()->GetXmin(),  hdists[ii]->GetXaxis()->GetXmax(), kFALSE, hnbx*10); ll->Add(hfit);
      hfit->SetLineColor(kRed);
      TString htit(Form("%s %d ", hdists[ii]->GetName(), kfail)); 
      for(Int_t ipar=0; ipar<npar; ipar++){
        htit+=Form("%e , ",pars[ipar]);
      }
      hfit->SetTitle(htit);
      hdists[ii]->SetTitle(htit);
      printf("%s\n", htit.Data());
      */
    }

    const Double_t fy = pars[1];
    const Double_t rmsy = TMath::Abs(pars[2]);
    
    const Double_t sysel= rmsy; //fy-ymin[ii];
    const Double_t syseh= rmsy; //ymax[ii]-fy;
    
    //they are already the truncated mean 0.1-0.9
    const Double_t statel = eyl[ii]; //eyl[ii]/nread;
    const Double_t stateh = eyh[ii]; //eyh[ii]/nread;
    
    printf("XGLUtils::GetGrSysMode0 setting %d %e %e %e %e -- %e %e -- gmin %d gmax %d\n", ii, xx[ii], fy, sysel, syseh, ymin[ii], ymax[ii], gmin[ii], gmax[ii]);
    
    grstat->SetPoint(ii, xx[ii], fy);
    grstat->SetPointError(ii, exl[ii], exh[ii], statel, stateh);

    grsys->SetPoint(ii, xx[ii], fy);
    grsys->SetPointError(ii, exl[ii], exh[ii], sysel, syseh);
  }
}


TH1D * XGLUtils::GetCDF(const TH1D *hh)
{
  const Int_t nbin = hh->GetNbinsX();
  const Double_t nor = hh->Integral(0,nbin+1);
  if(nor<EPSILON){
    printf("XGLUtils::GetCDF warning hh nor<epsilon! %e\n", nor);
    return 0x0;
  }

  TH1D *hcdf=(TH1D*)hh->Clone(Form("%scdf", hh->GetName()));
  for(Int_t ii=0; ii<=nbin+1; ii++){
    const Double_t cc=hh->Integral(0,ii);

    const Double_t eff = cc/nor;
    const Double_t err = TMath::Sqrt(eff*(1-eff)/nor);

    hcdf->SetBinContent(ii,eff);
    hcdf->SetBinError(ii,err);
  }

  hcdf->SetMaximum(1.01);

  return hcdf;
}

Int_t XGLUtils::GraphToData(const Int_t ng, TGraphAsymmErrors * grs[], Double_t xmin[], Double_t xmax[],  Double_t xs[], Double_t ys[], Double_t xes[], Double_t yes[])
{
  printf("\n");

  Int_t id=0;
  for(Int_t ig=0; ig<ng; ig++){
    printf("XGLUtils::GraphToData graph %d %s\n", ig, grs[ig]->GetName());
    const Int_t np = grs[ig]->GetN();

    for(Int_t ip=0; ip<np; ip++){
      Double_t xx=-999, yy=-999;
      grs[ig]->GetPoint(ip, xx, yy);

      if(xx<xmin[ig] || xx>xmax[ig])
        continue;

      Double_t ex = grs[ig]->GetErrorXlow(ip);
      Double_t ey = grs[ig]->GetErrorYlow(ip);
      xs[id] = xx;
      ys[id] = yy;
      xes[id] = ex;
      yes[id] = ey;
      printf("XGLUtils::GraphToData point %d %e %e %e %e\n", id, xx, yy, ex, ey);
      id++;
    }
  }
  printf("\n");

  return id;
}

TGraphAsymmErrors * XGLUtils::TreeToGraph(const TString tname, const TString cut, Char_t * sx, Char_t * sy, Char_t * sex0, Char_t * sex1, Char_t * sey0, Char_t * sey1, const Bool_t kXLogToLinear, const Bool_t kprint)
{
  TTree * rawtree=(TTree*)gDirectory->Get(tname);
  if(!rawtree){
    printf("XGLUtils::TreeToGraph no tree %s\n", tname.Data());exit(1);
  }

  TTree *cuttree = rawtree->CopyTree(cut);
  TTree * tmptree=(TTree*) cuttree->Clone(Form("tmptree%s%s",tname.Data(), sy));//necessary to avoid address confusion when executed repeatedly.

  Double_t vx = -1e10, vy = -1e10, vex0 =-1e10, vex1 = -1e10, vey0 = -1e10, vey1 = -1e10;
  tmptree->SetBranchAddress(sx, &vx);
  tmptree->SetBranchAddress(sy, &vy);
  if(strlen(sex0) && strlen(sex1)){
    tmptree->SetBranchAddress(sex0, &vex0);
    tmptree->SetBranchAddress(sex1, &vex1);
  }
  if(strlen(sey0)){
    tmptree->SetBranchAddress(sey0, &vey0);
  }
  //if sey0 and sey1 are the same, vey0 will not be set from tree, has to be set by vey1
  if(strlen(sey1)){
    tmptree->SetBranchAddress(sey1, &vey1);
  }
  
  TGraphAsymmErrors * tmpgr = new TGraphAsymmErrors;
  const Int_t ntot = tmptree->GetEntries();
  for(Int_t ii=0; ii<ntot; ii++){
    tmptree->GetEntry(ii);

    //remove trivial -1e10
    if(vex0<-1e9) vex0=0;
    if(vex1<-1e9) vex1=0;
    if(vey0<-1e9) vey0=0;
    if(vey1<-1e9) vey1=0;

    //minos error low < 0 by convention
    vex0=TMath::Abs(vex0);
    vex1=TMath::Abs(vex1);
    vey0=TMath::Abs(vey0);
    vey1=TMath::Abs(vey1);

    if(strcmp(sey0,sey1)==0){
      vey0=vey1;
    }

    if(kXLogToLinear){
      const Double_t xtrue = pow(10, vx);
      const Double_t xlow = pow(10, vx-vex0);
      const Double_t xhigh = pow(10, vx+vex1);
      vx = xtrue;
      vex0 = xtrue - xlow;
      vex1 = xhigh - xtrue;
    }

    tmpgr->SetPoint(ii, vx, vy);
    tmpgr->SetPointError(ii, vex0, vex1, vey0, vey1);
    if(kprint){
      printf("XGLUtils::TreeToGraph %s %s %s -- %d %e %e -- %e %e %e %e\n", tname.Data(), sx, sy, ii, vx, vy,  vex0, vex1, vey0, vey1);
    }
  }
  delete tmptree;

  TGraphAsymmErrors * gr = (TGraphAsymmErrors *)tmpgr->Clone(Form("gr%s%s",tname.Data(), sy));
  delete tmpgr;
  return gr;
}

TGraphAsymmErrors *XGLUtils::HistToGraph(const TString grname, const TH1 *hh, const Double_t thres, const TH1 *herr, const Double_t xmin, const Double_t xmax, const Bool_t klogx)
{
  const Int_t nbin = hh->GetNbinsX();
  Double_t xxs[nbin], yys[nbin], exls[nbin], exhs[nbin], eys[nbin];
  Int_t np=0;
  for(Int_t ii=1; ii<=nbin; ii++){
    const Double_t iyy = hh->GetBinContent(ii);
    if(iyy<=thres)
      continue;

    const Double_t iey = hh->GetBinError(ii);
    if(iey<EPSILON){
      if(iyy>EPSILON){
        printf("XGLUtils::HistToGraph warning! should be fine if this is ratio %d %e %e\n", ii, iyy, iey); //exit(1);
      }
      //continue;
    }

    const Double_t ixx = klogx? TMath::Sqrt(hh->GetXaxis()->GetBinLowEdge(ii)*hh->GetXaxis()->GetBinUpEdge(ii)): hh->GetBinCenter(ii);
    if(ixx<xmin || ixx>xmax){
      //printf("test XGLUtils::HistToGraph rejecting ixx %e xmin %e xmax %e\n", ixx, xmin, xmax);
      continue;
    }

    Double_t iexl = 0;
    Double_t iexh = 0;
    if(herr){
      iexl = iexh = herr->GetBinContent(herr->GetXaxis()->FindBin(ixx));
    }
    else{
      if(klogx){
        iexl = ixx - hh->GetXaxis()->GetBinLowEdge(ii);
        iexh = hh->GetXaxis()->GetBinUpEdge(ii) - ixx;
      }
      else{
        iexl = iexh = hh->GetBinWidth(ii)/2.;
      }

    }

    xxs[np] = ixx;
    yys[np] = iyy;
    exls[np] = iexl;
    exhs[np] = iexh;
    eys[np] = iey;
    np++;
  }
  TGraphAsymmErrors * gr = new TGraphAsymmErrors(np, xxs, yys, exls, exhs, eys, eys);
  gr->SetName(grname);
  gr->SetMaximum(hh->GetMaximum());
  gr->SetMinimum(hh->GetMinimum());
  gr->GetXaxis()->SetLimits(hh->GetXaxis()->GetXmin(), hh->GetXaxis()->GetXmax());
  return gr;
}

TH1D * XGLUtils::GetSeparation(const TH1D *hmpv, const TH1D *hwid, const Int_t anchor, const Bool_t kprint, const Double_t gr0)
{
  const Int_t nx = hmpv->GetNbinsX();
  const Double_t xmin = hmpv->GetXaxis()->GetXmin();
  const Double_t xmax = hmpv->GetXaxis()->GetXmax();

  TH1D *hsep = new TH1D(Form("%ssep",hmpv->GetName()), "", nx,  xmin, xmax);

  const Double_t m0 = hmpv->GetBinContent(anchor);
  const Double_t e0 = hmpv->GetBinError(anchor);
  const Double_t r0 = gr0>0 ? gr0 : hwid->GetBinContent(anchor);

  if(kprint) printf("XGLUtils::GetSeparation m0 %f ; r0 %f\n", m0,r0);

  for(Int_t ii=0; ii<=hsep->GetNbinsX()+1; ii++){
    hsep->SetBinContent(ii,0);
    hsep->SetBinError(ii,0);

    const Double_t mi = hmpv->GetBinContent(ii);
    const Double_t ei = hmpv->GetBinError(ii);
    const Double_t ri = hwid->GetBinContent(ii);
    if(ei<EPSILON){
      if(mi>EPSILON){
        printf("XGLUtils::GetSeparation error! %d %e %e\n", ii, ei, mi); exit(1);
      }
      continue;
    }

    const Double_t cont = fabs(mi-m0)/(r0+ri)*2;
    const Double_t err = sqrt(e0*e0+ei*ei)/(r0+ri)*2;
    if(kprint) printf("XGLUtils::GetSeparation %f: mi %f ; ri %f ==> %f\n",hmpv->GetBinCenter(ii), mi,ri, cont);
 
    hsep->SetBinContent(ii, cont);
    hsep->SetBinError(ii,err);
  }

  return hsep;
}

void XGLUtils::ToGraph(const TString fintxt, const Int_t istart, const Int_t nvar, TGraph * gr[])
{
  ifstream fr(fintxt);
  if(!fr){
    printf("XGLUtils::ToGraph error file not open!! %s\n", fintxt.Data()); exit(1);
  }

  const Int_t maxnpnt = 200;
  Double_t xx[maxnpnt];
  Double_t yy[nvar][maxnpnt];
  Int_t npnt = 0;

  Char_t buff[500];
  while(fr.getline(buff, 500)){
    TString pnum(buff);
    if(pnum.Contains("#"))
      continue;

    //---------- remove empty space in the beginning ---------
    while(pnum[0]==' '){
      pnum = pnum( 1, pnum.Length());
    }

    //---------- remove words in front of istart ------------
    for(Int_t ii=0; ii<istart; ii++){
      pnum = pnum( pnum.First(" ")+1, pnum.Length());
    }

    //---------- read x, remove 1 word ----------
    sscanf(pnum.Data(),"%lf",&(xx[npnt]));
    pnum = pnum( pnum.First(" ")+1, pnum.Length());

    //---------- for each y, read y remove 1 word ----------
    for(Int_t ivar=0; ivar<nvar; ivar++){
      sscanf(pnum.Data(), "%lf", &(yy[ivar][npnt]));
      pnum = pnum( pnum.First(" ")+1, pnum.Length());
    }

    npnt++;
    if(npnt>maxnpnt){
      printf("XGLUtils::ToGraph error npnt>maxnpnt!! %d %d\n", npnt, maxnpnt); exit(1);
    }
  }

  for(Int_t ivar=0; ivar<nvar; ivar++){
    //--------- Set Maximum of gr -----------
    Double_t ming = 1e10, maxg = -1e10;
    for(Int_t ipnt=0; ipnt<npnt; ipnt++){
      if(maxg<yy[ivar][ipnt])
        maxg = yy[ivar][ipnt];
      if(ming>yy[ivar][ipnt])
        ming = yy[ivar][ipnt];

      printf("XGLUtils::ToGraph Please Check %d %d: %e %ee\n", ivar, ipnt, xx[ipnt], yy[ivar][ipnt]);
    }
    printf("\n");

    gr[ivar] = new TGraph(npnt, xx, yy[ivar]);
    
    gr[ivar]->SetMinimum(ming);
    gr[ivar]->SetMaximum(maxg);
  }
}


TH1D* XGLUtils::ToPDF(const TH1 *hraw, const TString hn)
{
  TH1D * hist = (TH1D*) hraw->Clone((hn+hraw->GetName())+"pdf");
  Double_t ntot = 0;
  hist->Scale(0);

  //only consider region in range, ignore over/underflow. Important for fitting experiment data where under/overflow is ignored.
  const Int_t x0 = hist->GetXaxis()->GetFirst();
  const Int_t x1 = hist->GetXaxis()->GetLast();

  for(Int_t ii=x0; ii<=x1; ii++){
    const Double_t err = hraw->GetBinError(ii);
    const Double_t cont = hraw->GetBinContent(ii);

    //skip empty bins
    if(err<EPSILON){
      if(cont>EPSILON){
        printf("XGLUtils::ToPDF error! %d %e %e\n", ii, err, cont); exit(1);
      }
      continue;
    }

    const Double_t nn = cont*cont/err/err;
    hist->SetBinContent(ii,nn);
    hist->SetBinError(ii, TMath::Sqrt(nn));
    ntot += nn;
  }

  if(ntot<EPSILON){
    printf("XGLUtils::ToPDF ntot<epsilon ! %f %f\n", hist->GetEntries(), ntot);
    exit(1);
  }
  
  const Double_t bw = hist->GetBinWidth(1);
  for(Int_t ib=x0; ib<=x1; ib++){
    const Double_t cont = hist->GetBinContent(ib);
    if(cont<EPSILON)
      continue;

    //in case of finit number of bins (i.e. eff not always small), Binomial error is more accurate than Poisson error
    const Double_t eff = cont/ntot;
    const Double_t pdf = eff/bw;

    const Double_t dpdf = sqrt(eff*(1-eff)/ntot) / bw;
    hist->SetBinContent(ib, pdf);
    hist->SetBinError(ib, dpdf);
  }

  hist->SetEntries(ntot);

  return hist;
}

void XGLUtils::PrintOut(Int_t & id, const TString pn, const TCanvas *cc, const Int_t kopt)
{
  PrintOut(Form("ID%03d_%s", id++, pn.Data()), cc, kopt);
}

void XGLUtils::PrintOut(const TString pn, const TCanvas *cc, const Int_t kopt)
{
  if(!kopt)
    return;

  if(kopt & KPNG)
    cc->Print(pn+".png");
  if(kopt & KPDF)
    cc->Print(pn+".pdf");
  if(kopt & KROOT)
    cc->Print(pn+".root");
  if(kopt & KEPS)
    cc->Print(pn+".eps");
  if(kopt & KC)
    cc->Print(pn+".C");
  if(kopt & KJPG)
    cc->Print(pn+".jpg");
  if(kopt & KPS)
    cc->Print(pn+".ps");
}

TH2D* XGLUtils::NormalHist(const TH2D *hraw, const Double_t thres, const Bool_t kmax)
{
  //
  //normalized histogram
  //

  TH2D *hh=(TH2D*)hraw->Clone(Form("%snor",hraw->GetName()));
  hh->Scale(0);

  const Int_t x0 = hh->GetXaxis()->GetFirst();
  const Int_t x1 = hh->GetXaxis()->GetLast();
  const Int_t y0 = hh->GetYaxis()->GetFirst();
  const Int_t y1 = hh->GetYaxis()->GetLast();

  Double_t hmax = -1e10;
  Double_t hmin = 1e10;
  Double_t nent = 0;
  for(Int_t ix=x0; ix<=x1; ix++){

    //if option "e" is specified, the errors are computed. if option "o" original axis range of the taget axes will be kept, but only bins inside the selected range will be filled.
    TH1D * sliceh = hraw->ProjectionY(Form("tmpnormalhist%sx%d", hh->GetName(), ix), ix, ix, "oe");
    const Double_t tot = sliceh->GetEntries();

    TH1D * pdfh=0x0;

    if(tot>EPSILON){
      nent += tot;
    
      Double_t imax = -999;

      if(!kmax){
        pdfh = ToPDF(sliceh,"tmp");
      }
      else{
        imax = sliceh->GetBinContent(sliceh->GetMaximumBin());
      }

      for(Int_t iy=y0; iy<=y1; iy++){
        const Double_t cont = kmax ? sliceh->GetBinContent(iy)/imax : pdfh->GetBinContent(iy);
        const Double_t ierr = kmax ? sliceh->GetBinError(iy)/imax   : pdfh->GetBinError(iy);
        if(tot>thres && cont>0){
          hh->SetBinContent(ix, iy, cont);
          hh->SetBinError(ix,iy, ierr);
          if(cont>hmax) hmax = cont;
          if(cont<hmin) hmin = cont;
        }
      }
    }

    delete pdfh;
    delete sliceh;
  }

  hh->SetEntries(nent);
  hh->SetMinimum(0.99*hmin);
  hh->SetMaximum(1.1*hmax);
  return hh;
}

Double_t * XGLUtils::GetAxisArray(TAxis * aa)
{
  TArrayD * xs = new TArrayD(*(aa->GetXbins()));
  if(xs->GetSize()==0){
    const Int_t nbin=aa->GetNbins();
    xs->Set(nbin+1);
    for(Int_t ii=0; ii<=nbin; ii++){
      xs->SetAt(aa->GetBinUpEdge(ii),ii);
    }
  }

  Double_t * bins = xs->GetArray();
  if(!bins){
    printf("XGLUtils::GetAxisArray bins null!! %d\n", xs->GetSize());exit(1);
  }
  return bins;
}

void XGLUtils::FitSlicesY(const TH2D *hh, FFunc ffunc, TH1D *&hnor, TH1D *&hmpv, TH1D *&hwid, TH1D *&hres, TH1D *&hchi, const Double_t thres, const Double_t lowfrac, const Double_t highfrac, TList *ll)
{
  //
  //ffunc 3 parameter function: ALandau, AGaus
  //

  const Int_t x0 = hh->GetXaxis()->GetFirst();
  const Int_t x1 = hh->GetXaxis()->GetLast();
  const Int_t y0 = hh->GetYaxis()->GetFirst();
  const Int_t y1 = hh->GetYaxis()->GetLast();

  const Int_t nx = hh->GetNbinsX();  
  const Int_t ny = hh->GetNbinsY();
  const Double_t xmin = hh->GetXaxis()->GetXmin();
  const Double_t xmax = hh->GetXaxis()->GetXmax();
  const Double_t ymin = hh->GetYaxis()->GetXmin();
  const Double_t ymax = hh->GetYaxis()->GetXmax();

  hnor = new TH1D(Form("%s_amp",hh->GetName()), "", nx, xmin, xmax);
  hmpv = new TH1D(Form("%s_mpv",hh->GetName()), "", nx, xmin, xmax);
  hwid = new TH1D(Form("%s_wid",hh->GetName()), "", nx, xmin, xmax);
  hres = new TH1D(Form("%s_res",hh->GetName()), "", nx, xmin, xmax);
  hchi = new TH1D(Form("%s_chi",hh->GetName()), "", nx, xmin, xmax);

  const Double_t *hxbins = GetAxisArray(hh->GetXaxis());

  hnor->GetXaxis()->Set(nx, hxbins);
  hmpv->GetXaxis()->Set(nx, hxbins);
  hwid->GetXaxis()->Set(nx, hxbins);
  hres->GetXaxis()->Set(nx, hxbins);
  hchi->GetXaxis()->Set(nx, hxbins);

  for(Int_t ix=x0; ix<=x1; ix++){
    TH1D *htmp = new TH1D(Form("%s_%d", hh->GetName(), ix),"",ny, ymin, ymax);
    //checked, ok
    const Double_t *hhybins = GetAxisArray(hh->GetYaxis());
    //hh->GetYaxis()->GetXbins()->GetArray();
    if(hhybins){
      htmp->GetXaxis()->Set(ny, hhybins);
    }

    Double_t ntot = 0;
    for(Int_t iy=y0; iy<=y1; iy++){
      const Double_t be = hh->GetBinError(ix,iy);
      const Double_t bc = hh->GetBinContent(ix, iy);

      if(be<EPSILON){
        if(bc>EPSILON){
          printf("XGLUtils::FitSlicesY error %d %d %e %e\n", ix, iy, be, bc); exit(1);
        }
        continue;
      }

      htmp->SetBinContent(iy, bc);
      htmp->SetBinError(iy, be);

      ntot += (bc/be)*(bc/be);

      //if(be) printf("test %d %d : %f %f %f\n", ix, iy, bc, be, pow(bc/be,2));
    }

    hnor->SetBinContent(ix, ntot);
    hnor->SetBinError(  ix, 0);
    
    if(ntot<thres || htmp->GetRMS()<EPSILON){
      delete htmp;
      continue;
    }

    //test htmp->Draw();
    Double_t pars[10]={htmp->Integral(0,htmp->GetNbinsX()+1)*htmp->GetBinWidth(1), htmp->GetMean(), htmp->GetRMS()};
    Double_t errs[10], chi[10];
    
    TH1D * tmpslicefit = 0x0;
    if(ffunc){
      fgFitPrintLevel = -1;

      const Int_t kfail = ChisquareFit(htmp, ffunc, 3, pars, errs, chi);
      if(kfail){
        printf("XGLUtils::FitSlicesY ChisquareFit fail! %s %d n %.0f, direct continue\n", htmp->GetName(), kfail, htmp->Integral(0,htmp->GetNbinsX()+1));
        delete htmp;
        continue;
      }
      tmpslicefit = GetHfit(Form("%sfit", htmp->GetName()), ffunc, pars, htmp->GetXaxis()->GetXmin(), htmp->GetXaxis()->GetXmax(), kFALSE);
    }
    else{
      pars[1] = TruncatedMean(htmp, lowfrac, highfrac, &(pars[2]), &(errs[1]));
      errs[2]=0;
      chi[0] = chi[1] = 0;
    }

    pars[2]=TMath::Abs(pars[2]);
    //hnor->SetBinContent(ix, htmp->GetBinContent(htmp->GetMaximumBin()));//htmp->Integral(0,htmp->GetNbinsX()+1));
    hmpv->SetBinContent(ix, pars[1]);
    hmpv->SetBinError(  ix, errs[1]);

    hwid->SetBinContent(ix, pars[2]);
    hwid->SetBinError(  ix, errs[2]);

    hres->SetBinContent(ix, fabs(pars[1])>EPSILON? pars[2]/fabs(pars[1]):0);
    hres->SetBinError(  ix, fabs(pars[1])>EPSILON? errs[2]/fabs(pars[1]):0);

    hchi->SetBinContent(ix, chi[1]>=1 ? chi[0]/chi[1]: 0);
    hchi->SetBinError(ix, 0);

    if(ll){
      ll->Add(htmp);
      if(tmpslicefit) 
        ll->Add(tmpslicefit);
    }
    else{
      delete htmp;
      delete tmpslicefit;
    }
  }

  TH1 *hhs[]={hnor, hmpv, hwid, hres, hchi};
  const TString yt[]={"N", "MPV", "#sigma", "#sigma/MPV", "#chi^{2}/NDOF"};
  const Int_t nh = sizeof(hhs)/sizeof(TH1*);
  for(Int_t ii=0; ii<nh; ii++){
    hhs[ii]->SetYTitle(Form("%s of %s", yt[ii].Data(), hh->GetYaxis()->GetTitle()));
    hhs[ii]->SetXTitle(hh->GetXaxis()->GetTitle());
    hhs[ii]->GetYaxis()->SetTitleOffset(hh->GetYaxis()->GetTitleOffset());
    hhs[ii]->SetTitle(hh->GetTitle());
  }
}

Double_t XGLUtils::TruncatedMean(const TH1 *hh, const Double_t lowfrac, const Double_t highfrac, Double_t * grms, Double_t * gerr)
{
  //
  //do truncation on histogram
  //
  //if hh is scaled, be sure Sumw2 is called before scaling!! then mean, rms and err will all be correct
  
  //with under-/over-flow
  Double_t npreTrunc = 0;
  for(Int_t itmp=0; itmp<=hh->GetNbinsX()+1; itmp++){
    const Double_t be = hh->GetBinError(itmp);
    const Double_t bc = hh->GetBinContent(itmp);
    if(be<EPSILON){
      if(bc>EPSILON){
        printf("XGLUtils::TruncatedMean (hist) error %e %e %d\n", bc, be, itmp); exit(1);
      }
      continue;
    }
    npreTrunc += bc*bc/be/be;
  }

  const Double_t nstart = npreTrunc*lowfrac;
  const Double_t nstop = npreTrunc*highfrac;

  //with Double_t this should also handle normalized hist
  Double_t ntot = 0;
  Double_t nused = 0;
  Double_t sum = 0;
  Double_t w2s = 0;
  for(Int_t itmp=0; itmp<=hh->GetNbinsX()+1; itmp++){
    const Double_t be = hh->GetBinError(itmp);
    const Double_t bc = hh->GetBinContent(itmp);
    if(be<EPSILON){
      if(bc>EPSILON){
        printf("XGLUtils::TruncatedMean (hist) error %e %e %d\n", bc, be, itmp); exit(1);
      }
      continue;
    }
    const Double_t weight = bc*bc/be/be;
    ntot+=weight;
    //<= correct, because when high=1, nstop has to be included
    if(ntot>nstart && ntot<=nstop){

      const Double_t val = hh->GetBinCenter(itmp);
      sum += weight*val;
      w2s += weight*val*val;
    
      nused += weight;

      //printf("test %d %f %f --- %f %f -- %f %f\n", itmp, weight, val, sum, w2s, nused, nsample);
    }
    else if(ntot>nstop){
      if(itmp>=hh->GetNbinsX()){
        printf("XGLUtils::TruncatedMean warning hist range too small %s %f %f %d %d, %15f %15f %15f; nused w2s sum set to 0\n", hh->GetName(), hh->GetBinLowEdge(1), hh->GetBinLowEdge(itmp), itmp, hh->GetNbinsX(), hh->GetBinContent(hh->GetNbinsX())/hh->Integral(0,hh->GetNbinsX()+1), hh->GetBinContent(hh->GetNbinsX()), hh->Integral(0,hh->GetNbinsX()+1)); //exit(1);
        nused = 0;
        w2s = sum = 0;
      }
      break;
    }
  }

  return GetMeanRMS(nused, sum, w2s, grms, gerr);

}

//===================================================================
//                           IO
//===================================================================

TChain * XGLUtils::InputFiles(const TString file, const TString tr, Char_t *dir)
{
  TChain *ch=new TChain(tr);
  
  if(file.Contains(".root"))
    ch->Add(file);
  else{
    ifstream fin(file);
    if(!fin){
      printf("XGLUtils::InputFiles file not found \n%s\n\n",file.Data()); exit(1);
    }
    
    TString buff;
    while(fin.good()){
      fin>>buff;
      if(buff!=""){
        if(dir){
          buff.Prepend(dir);
        }
        ch->Add(buff);
      }
    }
  }

  //const Int_t ent = ch->GetEntries(); //takes infinity time!!
  printf("\t%d trees!\n",ch->GetNtrees());
  
  return ch;
}

//===================================================================
//                            Physics
//===================================================================

Double_t XGLUtils::ToPt(const Double_t mom, const Double_t eta)
{
  /*
    TLorentzVector b
    b.SetPtEtaPhiE(0.567,0.123,0,0)
    b.Pt()-XGLUtils::ToPt(b.P(), b.Eta())
    //should be eps
   */

  const Double_t sintheta = TMath::Sin( 2*TMath::ATan( TMath::Exp( TMath::Abs( eta ) ) ) );
  return mom*sintheta;
}

Double_t XGLUtils::Eta2Y(const Double_t mass, const Double_t pt, const Double_t eta)
{
  //checked!
  /*
    TLorentzVector b
    b.SetPtEtaPhiM(1, 0.7, 0, 1)
    double eta = 0.7, mass = 1, pt = 1
    const Double_t sinh = TMath::SinH(eta)/TMath::Sqrt(1+TMath::Power(mass/pt,2));
    TMath::ASinH(sinh) - b.Rapidity()
    b.SetPtEtaPhiM(1, 0.7, 3, 1)
    TMath::ASinH(sinh) - b.Rapidity()
    //should be eps
  */

  const Double_t sinhy = TMath::SinH(eta)/TMath::Sqrt(1+TMath::Power(mass/pt,2));
  return TMath::ASinH(sinhy);
}

Double_t XGLUtils::DeltaY(const Double_t mass, const Double_t pt, const Double_t eta0, const Double_t eta1)
{
  const Double_t y0 = Eta2Y(mass, pt, eta0);
  const Double_t y1 = Eta2Y(mass, pt, eta1);
  return TMath::Abs(y0-y1);
}

